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1.  INTRODUCTION 


In  a  recent  report  (1]  a  mixed  numericai/analytical  approach  to  the  computation  of  a  ring- 
symmetric  spacecraft  exhaust  plume  was  presented.  The  numerical  scheme  had  been  implemented  in 
a  code  named  "JET"  which  is  capable  of  generating  whole-plume  flow  fields,  while  the  analytic 
approximation  is  restricted  to  the  ring- symmetric  centered  rarefaction  waves  (CRW)  that  flank  the 
plume.  The  present  report  is  intended  to  serve  as  a  supplement  to  [1]  in  providing  details  on  the 
computational  scheme  and  the  code  JET. 

The  spacecraft  exhaust  flow  (Fig.  1  of  (11  )  is  idealized  as  a  ring- symmetric  steady  isentropic 
expansion  of  an  ideal  gas.  The  nozzle  lips  are  assumed  sharp;  the  supersonic  flow  from  the  exit 
surface  of  the  ring-nozzle  is  assumed  uniform,  and  the  background  is  considered  to  be  perfect 
vacuum. 

The  standard  scheme  for  computing  such  idealized  ring-plumes  is  the  classical  (direct)  method  of 
characteristics  [2] .  At  a  preliminary  phase  of  the  present  laser  exhaust  study,  a  code  AXS  YM  [3]  was 
written  for  computing  ring-plumes  using  this  method.  A  notorious  shortcoming  of  the  direct  method 
of  characteristics  is  that  the  solution  grid  is  highly  irregular,  being  formed  by  the  (oblique) 
intersection  of  the  C  +  and  the  C  ”  families  of  characteristic  lines.  We  first  encountered  a  difficulty 
with  this  grid  while  seeking  a  scheme  for  integrating  the  molecular  opacity  along  a  straight  line  [1]  . 
This  computation  would  have  required  rather  complex  coding  for  the  geometry  of  intersection 
between  a  straight  line  and  an  irregular  grid.  It  seemed  preferable  to  opt  for  a  computation  scheme 
that  would  produce  a  more  regular  grid,  even  at  the  expense  of  some  loss  of  accuracy.  Such  scheme 
is  the  inverse  marching  method  of  characteristics  [4]  . 

Generally  the  marching  in  this  type  of  scheme  is  in  the  downstream  direction,  i.e..  the  y  direction 
in  our  case.  Grid  points  are  located  on  a  succession  of  constant-y  rows,  thereby  introducing  a 
measure  of  regularity  in  the  solution  grid.  Just  two  rows  have  to  be  stored  in  the  computer  core 
memory  -  the  "old"  row  and  the  "new"  row,  whereas  in  the  direct  method  of  characteristics  whole 
grid-image  matrices  are  required  to  reside  simultaneously  in  core  memory'. 

The  first  version  of  the  JET  code  was  based  on  the  inverse  marching  scheme  given  by  Zucrow  and 
Hoffman  (Section  12-5  in  [4j  ),  where  the  flow  variables  were  the  two  cartesian  velocity  components.  \ 
The  computation  seemed  accurate  everywhere,  except  within  the  centered  rarefaction  wave  (CRW).  i 
In  an  attempt  to  replicate  a  planar  CRW  (Prandtl-Meyer  flow),  the  numerical  solution  exhibited  an 


instability  :  Mach  number  increased  along  the  (low  pressure)  boundary  characteristic  line,  rather  than 
remain  constant. 

A  qualitative  explanation  for  this  instability  is  the  following.  Flow  gradients  in  a  CRW  are 
inversely  proportional  to  distance  from  the  comer,  so  that  the  inverse  marching  scheme  gives  rise  to 
an  amplification  of  interpolation  errors  at  every  marching  step,  leading  to  an  apparently  divergent 
(unstable)  numerical  solution.  Increasing  the  order  of  interpolation  from  linear  to  cubic  did  not 
eliminate  the  instability. 

Looking  for  a  scheme  that  would  replicate  a  planar  CRW  accurately,  we  tried  the  modified 
marching  idea  as  presented  by  Zucrow  and  Hoffman  for  1-D  time-dependent  flows  (Sections  19-6(a) 
and  19-6(j)  of  (4] ).  In  this  scheme  new  grid  points  are  determined  by  forwardly  extending  a  "primary" 
family  of  continuous  characteristic  lines  from  old  grid  points.  The  primary  family  in  a  CRW  is  the 
characteristics  fanning  out  from  the  corner  (we  assume  it  is  the  C+  family).  By  choosing  this 
modified  scheme,  the  interpolation  for  trace  points  obtained  from  reversely  extended  C  +  lines  was 
eliminated.  However,  the  corresponding  interpolation  for  the  transverse  C —  characteristics 
remained,  and  with  it  the  aforementioned  instability. 

In  order  to  replicate  a  planar  CRW,  we  had  to  replace  the  flow  variables  by  the  Riemann 
invariants  (v±0).  In  a  C+  planar  CRW,  the  Riemann  invariant  (v  +  0)  is  uniformly  constant,  so 
that  the  interpolation  in  (v  +  0)  due  to  reversely  extending  C  -  characteristics  introduces  no  error  at 
all.  This  scheme,  which  we  named  SIMA  (Semi  Inverse  Marching  Algorithm),  was  indeed  verified  to 
replicate  a  planar  CRW  exactly,  when  implemented  in  the  code  JET. 

The  plan  of  this  report  is  the  following.  In  Ch.  2  w’e  supplement  the  description  of  the  numerical 
scheme  given  in  Ch.  2  of  [1]  ,  by  adding  more  details  on  the  computational  procedure.  A  description 
of  the  code  JET  is  given  in  Ch.  3,  and  the  code  listing  is  reproduced  in  Ch.  4. 

Note  on  symmetry  : 

The  code  JET  has  two  symmetry  options.  When  DELTA  =1  a  ring-symmetry  is  in  effect:  when 
DELTA  =  0.  a  planar  symmetry  is  in  effect.  An  axisymmetric  jet  exiting  in  the  y  direction  from  the 
same  nozzle  aperture  along  the  x  axis  can  readily  be  computed  by  replacing  all  terms  in  the  code  that 
correspond  to  sin(0)/y  in  the  compatibility  equations  (2.1-1),  by  cos (0)/x.  In  that  case  the  coding  is 
virtually  unchanged,  and  the  only  care  that  should  be  exercised  is  for  the  difference  equations  for  new 
grid  points  on  or  near  the  y  axis.  Also,  all  reference  to  the  analytic  approximation  of  the  ring- 
symmetric  CRW  [1)  should  be  deleted  in  this  case,  as  it  is  designed  specifically  for  ring-symmetry. 


2.  THE  COMPUTATIONAL  SCHEME 

A  basic  description  of  the  semi  inverse  (SIMA)  and  inverse  marching  schemes  was  given  in  Ch.  2 
of  [1]  .  We  supplement  this  description  by  specifying  the  slightly  modified  definition  of  Riemann 
invariants  in  the  code,  and  by  giving  information  about  some  ancillary  computations. 

2.1  Riemann  Invariants 

The  compatibility  equations  whose  integration  constitutes  the  numerical  solution  to  the  governing 
equations  [1]  are  expressed  in  terms  of  the  Riemann  invariants  as  follows : 

Along  C  +  -  (v  -  0)4  -  (v  -  0)2  +  to  sinji24  sin024  An  /  y24 

(2.1-1) 

Along  C-  -  (v  +  0)4  -  (v  +  0)j  +  to  sinfi,4  sin014  A^  /  y14 

The  Riemann  invariants  (v  ±  0)  are  modified  for  convenience,  by  adding  a  constant  to  both  v  and 
0.  The  new  definitions  of  v(M)  and  0  are  : 

v(M)  =  -T  arctan(Tq)+  arctan(q) 

q  =  (M2-l)‘l/2  (2.1-2) 

0  — ♦  0  ~  0^ 

Thus,  in  a  Prandtl-Meyer  flow  with  entry  Mach  number  of  M,,  the  modified  values  of  both  v(M) 
and  0  vanish  as  M-»ao.  As  a  consequence,  in  a  C+  Prandtl-Meyer  flow  the  modified  invariant 
(v  +  0)  vanishes  uniformly.  In  this  modified  form,  the  computation  of  M  from  v(M)  is  readily  done 
by  performing  standard  Newton- Raphson  iterations  (in  RFL'NC),  using  the  derivative  : 

v'(q)  =  -(r-uIO  +  rVxi+q2)]'1  (2.1-3) 
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2.2  The  Integration  Scheme  for  a  New  Grid  Point 


£ 


The  integration  scheme  has  been  sketched  in  Ch.  2  of  [1] .  It  is  performed  in  INVMAR  for  inverse 
marching  points  or  in  SEMINV  for  semi-inverse  marching  (SIMA)  points.  The  computational 
scheme  is  specified  via  the  following  seven-step  procedure  : 


INVMAR  (Inverse  Marching) 

(a)  Grid  :  At  this  stage  the  new  grid  point  has  already  been  defined. 

(b)  Predictor  :  Flow  variables  are  the  interpolated  (linear  nearest-neighbor)  value  on  the  old  row 
for  a  point  having  the  new  grid  x  coordinate  (x4). 

(c)  Centered  variables  :  Denote  the  Riemann  invariants  by 


RM  =  (v  +  0) 
RP  =  (v  — 0) 


(2.2-1) 


then  centered  values  for  segments  (1,4)  and  (2,4)  (using  code  notation)  are  : 


RM14  =  (RMl  +  RM4)/2 
RM24  =  (RM2  +  RM4)/2 


RP14  =  (RP1  +  RP4)/2 
RP24  *  (RP2  +  RP4)/2 


(2.2-2) 


All  other  centered  flow  variables  are  computed  from  the  centered  Riemann  invariants  by  calling 
RFL'NC. 

(d)  Inverse  Extension  :  old  trace  points  x.,  x.,  are  evaluated  from  the  geometrical  relations 
Along  C-  ....  Vnew  -  yold  =  (X4-Xj)  tan(0u-nu) 

(2.2-3) 

Along  C  +  ....  ynew  -  yo|d  =  (x4  -  X2)  tan(024  +  p,4) 


(e)  Interpolation  :  find  Riemann  invariants  RM,  RP  at  old  trace  points  x,  and  X,  through  nearest- 
neighbor  linear  interpolation  by  calling  INTERP. 

(0  Integration  :  Using  the  compatibility  relations  in  finite-difference  form  (2.1-1)  with  segment- 
centered  coefficients,  compute  iteration-updated  values  of  Riemann  invariants  at  new  grid  point. 


(g)  Corrector  :  if  values  of  Riemann  invariants  and  old  trace  points  Xj,  are  not  sufficiently 
convergent,  resume  the  procedure  at  step  (c)  above. 

SEMINV  (Semi  Inverse  Marching  -  SIMA) 

(a)  Grid  :  New  grid  point  (x4)  is  determined  as  part  of  the  SIMA  scheme  at  step  (d)  below. 

(b)  Predictor  :  Flow  variables  are  those  of  point  (Xj.y^j). 

(c)  Centered  variables :  Identical  to  step  (c)  above. 

(d)  Semi-Inverse  Extension  :  new  grid  point  x4  and  old  trace  point  x,  are  evaluated  from  the 
geometrical  relations  in  Eq.  (2.2-3)  above. 

(e)  Interpolation :  find  Riemann  invariants  RM,  RP  at  old  trace  point  x,  through  nearest- 
neighbor  linear  interpolation  by  calling  INTERP. 

(0  Integration  :  Identical  to  step  (f)  above. 

(g)  Corrector  :  Identical  to  step  (g)  above,  except  for  replacing  x^  in  the  convergence  test  by  x4. 

2.3  Boundary  Conditions 

On  the  vacuum  side  the  boundary  conditions  (p  =  0)  can  only  be  approximately  implemented  in  a 
method  of  characteristics  scheme.  We  do  so  by  ending  the  computation  on  a  certain  ''final"  C  +  fan 
characteristic  line  that  starts  out  with  a  sufficiently  high  Mach  number  Mf  at  the  corner  (typically 
Mf=  34).  The  marching  computation  of  new  grid  points  on  the  boundary  C  characteristic  via  the 
SI  VIA  scheme  is  identical  to  that  of  C+  characteristics  within  the  ring- symmetric  CRW.  It  is  noted 
that  under  this  boundary  scheme  some  outflow'  takes  place  through  the  boundary  characteristic  line, 
so  that  the  total  mass  flow  through  a  row  y=  y„„„,  decreases  slightlv  as  v„.„  increases. 

At  the  nozzle  exit  the  boundarv  conditions  are  assumed  to  be  uniform  outflow  in  the  radial  (v) 
direction  with  Mach  number  Mj.  At  the  nozzle  lip,  the  SIMA  integration  starts  out  from  a  presumed 
planar  CRW  (Prandtl- Meyer  flow)  at  the  corner  (i.e.,  the  associate  CRW  in  the  terminology  of  Ch.  3 
in  [1)  ). 

At  the  plane  of  symmetry  (X  =  0)  the  boundary  condition  is  simply  Q  =  n/2.  However,  this 
condition  is  implemented  indirectly,  by  assuming  that  the  flow  at  virtual  grid  points  with  x  <  0  is  a 
mirror-image  of  the  flow  at  the  corresponding  X  >  0  points.  The  reason  is  that  when  a  new  grid 
point  of  X4  =  0  or  of  x4  sufficiently  close  to  zero  is  considered  for  inverse-marching  integration,  the 
inversely  extended  trace  point  (Xj ,yojd)  can  be  at  x  <  0.  Considering  the  subtraction  of  0L  from  0  as 
in  Eq.(2.1-2),  the  reflection  rules  are  : 


RM  -»  RP  +  (re  —  20l) 


(2.3-1) 


Cv: 


R. 


A 


RP 


RM  -  (n  -  20l) 


where  values  on  the  left  and  right  of  the  -♦  symbol  correspond  to  values  left  and  right  of  x  =  0.  This 
boundary  condition  is  implemented  in  INTERP. 


2.4  Continuum  Breakdown  Surface 


As  an  informative  option,  the  code  JET  can  compute  (in  PLUMES)  points  on  a  surface  of 
continuum  breakdown  (5,6,71 ,  which  is  defined  as  a  line  of  constant  B,  where  B  is  given  by  : 


B  =  -  (u/9)  p*‘  (dp/dS) 


-1/2 

q>  -  4{jry)  <r  n  a 


(2.4-1) 


When  the  standard  isentropic  relations  for  p  and  n  in  terms  of  M  are  substituted  in  (2.4-1),  the 
flow  speed  is  expressed  as  U  =  Ma  and  the  streamwise  gradient  of  M  is  expressed  in  cartesian 
coordinates,  we  get  : 


B  =  X0  (ny/8)l/2  M2  [l  +  ((y-l)/2)M‘] 


2i1/(Y-D-I 


[Mxcos0  +Mysin0j 


=  (2 1/2  <T  n0)  ‘ 


(2.4-2) 


Note  that  the  sign  of  B  has  been  chosen  as  positive  for  expansion  flows.  This  definition  is 
preferred  to  taking  an  absolute  value  of  the  flow  gradient,  since  it  assures  proper  interpolation  of  B 
even  if  its  spatial  distribution  goes  through  B  =  0. 


Due  to  the  dependence  of  B  on  a  spatial  gradient,  its  numerical  evaluation  (see  BREAK)  is 
attributed  to  mid-grid  points  both  in  x  and  in  y. 


3.  THE  JET  CODE 


In  this  chapter  we  provide  a  concise  description  of  the  JET  code  according  to  its  version  at  the 
time  of  the  JET018  run.  This  description  is  intended  as  an  aid  in  reading  the  code  listing  which  is 
given  in  Ch.  4. 

The  plan  of  this  chapter  is  as  follows.  Array  variables  that  constitute  the  mainstay  of  the 
computational  scheme  are  described  in  Section  3.1.  Auxiliary  array  variables  that  are  used  primarily 
for  processing  the  information  generated  by  the  numerical  scheme,  are  described  in  Section  3.2, 
followed  in  Section  3.3  by  a  list  of  major  parameters  that  control  the  computation  (some  of  them  also 
serve  as  run  data).  Finally,  all  subroutines  are  listed  and  described  in  Section  3.4. 


3.1  Main  Variables 

The  array  variables  used  for  the  computational  scheme  are  organized  in  two  labeled  COMMON 
groups.  The  first  group  /VECS /  is  designed  to  hold  two  grid  rows  -  the  old  row  designated  by  suffix 
F  and  the  new  row  designated  by  suffix  N.  The  second  group  /CHARAC/  are  characteristic-indexed 
arrays  that  hold  information  about  continuous  characteristic  lines.  This  characteristic  information  is 
used  in  two  ways  :  it  is  incorporated  in  the  SIMA  computational  scheme  for  the  CRW  region,  and  it 
is  used  to  store  data  for  optional  plotting  of  characteristic  lines  (see  PLLMES  and  PRINT). 


The  basic  organization  is  that  the  new  arrays  (suffix  N)  are  those  in  which  values  are  stored  during 
the  course  of  the  marching  computational  procedure.  At  the  end  of  each  marching  step,  values  are 
transferred  from  new  arrays  to  old  arrays  (suffix  F);  this  is  done  in  MOVE.  In  the  array  listing 
below,  we  indicate  in  parenthesis  the  subroutine  (or  subroutines)  in  which  that  new  array  is  defined. 


/VECS  / 

XN(I)  x  coordinate  of  grid  point  I.  (GRIDN) 

RMN(I)  modified  Riemann  invariant  (v  +  0)  at  grid  point  I.  (BEGIN.  INVMAR. 

LOADC). 

RPN(I)  modified  Riemann  invariant  (v-0)  at  grid  point  I.  (BEGIN,  INVMAR. 

LOADC). 


MN(I) 

ML'N(I) 

TETAN(I) 


Mach  number  at  grid  point  I  (BEGIN,  INVMAR,  LOADC). 

Mach  angle  p  at  grid  point  I.  (BEGIN,  INVMAR.  LOADC). 

true  (unmodified)  flow  angle  0  at  grid  point  I.  (BEGIN.  INVMAR,  LOADC). 


MUMMamra 


BN(I)  value  of  breakdown  parameter  B  at  point  1-1/2  (and  at  half  a  marching  step 

back  in  y  as  well).  (BREAK). 

XTEMP(I)  used  for  auxiliary  computation  of  1-1/2  grid  points  in  PLUMES. 


/CHARAC/ 

XCHARN(KC) 

YCHARN(KC) 

RMCARN(KC) 

RPCARN(KC) 

TCHARN(KC) 

MUCARN(KC) 

CSIGNN(KC) 

MCHARX(KC) 

MCHARI(KC) 


x  coordinate  of  point  on  characteristic  line  number  KC.  (BEGIXr, 
SEMIXV,  PLUMES). 

y  coordinate  of  point  on  characteristic  line  number  KC.  (BEGIN,  SEMIXV, 
PLUMES). 

modified  Riemann  invariant  (v  +  0)  of  point  on  characteristic  line  number 
KC.  (BEGIN,  SEMIXV). 

modified  Riemann  invariant  (v-0)  of  point  on  characteristic  line  number 
KC.  (BEGIN,  SEMIXV). 

true  (unmodified)  flow  angle  0  at  point  on  characteristic  line  number  KC. 
(BEGIN,  SEMIXV). 

Mach  angle  ji  at  point  on  characteristic  line  number  KC.  (BEGIN, 
SEMIXV). 

sign  of  characteristic  line  number  KC.  It  has  value  1  for  C  +  and  value  -  1 
for  C~ .  Note  that  upon  reflection  of  a  C  +  line  from  the  symmetry  plane 
(X=  0),  the  sign  value  is  changed  from  1  to  -  1.  (BEGIN,  SEMIXV). 

Mach  number  at  point  on  characteristic  line  number  KC.  (BEGIN. 
SEMIXV). 

Mach  number  at  Prandtl-Meyer's  fan  characteristic  number  KC  at  the 
comer.  It  is  defined  initially  and  is  not  changed  during  the  run.  (BEGIN). 


3.2  Auxiliary  Variables 

In  addition  to  the  major  arrays  mentioned  above,  there  are  several  groups  of  auxiliary’  arrays  that 
do  not  affect  the  computational  scheme,  but  are  intended  for  informative  processing  of  the  results. 
These  groups  are  /PLUME/,  /IPLUME /,  /THICKY/.  /THICKX/.  /GRP/.  /PLUME/  is  used  to 
preserve  points  on  special  lines  for  later  plotting  (in  a  separate  code).  /THICKY/  and  /THICKX/ 
are  for  storing  values  of  radial  { y  >  and  lateral  (x)  molecular  opacities.  The  group  /GRP/  is  used  in 
conjunction  with  comparative  computation  of  the  ring- symmetric  CRW  flow  according  to  the  analytic 
approximation  (1}  . 


yywkj v.uv 
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/PLUME/  (PLUMES,  PRINT) 

XPL(J,IPL)  X  coordinate  at  marching  step  J  of  special  line  number  IPL. 

YPL(J,IPL)  y  coordinate  at  marching  step  J  of  special  line  number  IPL. 

/IPLUME/  (PLUMES,  PRINT) 

KPL  number  of  special  lines  computed  in  PLUMES. 

ITYPL(IPL)  index  indicating  the  type  of  special  line  number  IPL. 

/THICKY/  (OPACY,  PRINT) 

XTH(J)  X  coordinate  on  boundary  characteristic  line  at  marching  step  J,  from  which 

radial  opacity  is  integrated. 

TH(J)  radial  opacity  computed  by  y-integration  from  the  boundary  point  defined  by 

XTH(J)  (up  to  current  YN). 

/THICKX/  (OPACX,  PLUMES,  PRINT) 

YXI(JXI)  y  coordinate  of  printed  row  number  JXI  (the  index  JXI  counts  just  rows  that 

have  been  printed).  The  row  to  be  printed  next  upon  calling  PRINT  is  the  row 
having  YF  near  YXI(JXI). 

XI(I,JXI)  lateral  (x)  molecular  opacity  [1]  at  point  XF(I),  for  printed  row  JXI.  It  is 

obtained  by  numerically  integrating  the  solution  obtained  from  the  JET 
computation  (see  OPACX). 

XIPM(I,JXI)  same  as  XI(I,JXI)  except  that  the  Prandtl-Meyer  solution  is  used  to  estimate 
the  flow  at  grid  points  XF(I). 

XIGRP(IJXI)  same  as  above,  except  that  the  analytic  approximation  to  a  ring-symmetric 
CRW  [l]  is  used  to  estimate  the  flow  at  grid  points  XF(I). 

XIAPP(I,JXI)  same  as  XIGRP(I,JXI)  except  that  the  numerical  integration  is  replaced  by  an 
approximate  closed-form  expression  [1]  . 

XIF(IJXI)  stores  grid  points  XF(1)  of  printed  row  JXI. 


/GRP /  (PRINT,  HMSET,  MFUNC,  HINTER,  MATCH) 


DMINV 

MHINV(I) 


HMV(I) 


increment  of  inverse  Mach  number  for  array  MHINV(I). 

inverse  Mach  number  array  (from  0  to  1/M  EXIT),  from  which  the  H(M) 

function  can  be  evaluated  (HMSET). 

values  of  the  H(M)  function  evaluated  by  numerical  integration.  It  is  used  to 
compute  this  function  by  interpolation.  (HMSET,  HINTER). 


3.3  Major  Parameters 


Parameters  that  define  and  control  a  particular  run  (such  as  the  maximum  y  for  the  marching 
computation,  the  number  of  grid  points  on  a  row  and  many  more)  are  defined  in  1NIDAT.  (The 
code  JET  has  no  input  file  and  no  READ  statements).  The  major  control  parameters  are  grouped  in 
/PAR/  (floating  point)  and  in  /IPAR/  (integers);  thermodynamic  data  are  grouped  in  /STAG/. 

We  indicate  in  the  listing  the  subroutines  in  which  the  labeled  COMMON  group  or  a  particular 
parameter  is  defined  (or  sometimes  referred  to). 

/PAR/  (INIDAT) 

MEXIT  nozzle  exit  Mach  number  (M,). 

MFIN  Mach  number  of  the  final  (boundary)  CRW  characteristic  at  the  comer  (Mf). 

YMAX  maximum  value  of  y  for  the  marching  scheme.  When  YF.GE.YMAX  the  run  is 
terminated. 

DYO  initial  marching  step. 

DY  current  marching  step. 

DYNEXT  next  marching  step  (YSTEP). 

STAB  stability  coefficient  for  marching  step  (STAB.LE.l).  (See  YSTEP). 

DELTA  symmetry  index.  DELTA  =  0  for  plane  symmetry;  DELTA*  1  for  ring-symmetry. 

PSIl  angle  of  Prandtl-Meyer  fan  characteristic  at  exit  conditions  (measured  from  x 

axis). 

PSIF  angle  of  final  (boundary)  Prandtl-Meyer  fan  characteristic. 

SIGMA  collision  cross-section  (<x). 

FRACG  the  number  of  intervals  initially  allocated  to  the  CRW  fan  is  a  FRACG  fraction  of 
the  total  number  of  intervals  (KFO-I).  (see  BEGIN). 

EPSIL  convergence  parameter  (small  number).  (IN'VMAR,  SEMIN'V). 

TETLIYI  flow  angle  (from  X  axis)  of  the  limiting  (p  =  0)  velocity  vector  of  the  flow  at  the  lip- 
centered  Prandtl-Meyer  fan. 

TETSYM  PAI-2*TETLIM  for  reflection  transformation  (see  INTERP). 

/IPAR/  (INIDAT) 

JMAX  maximum  number  of  marching  steps.  If  J.GE.JMAX  run  is  terminated. 

KFO  initial  (and  maximum)  number  of  grid  points  in  a  row. 

KF  current  number  of  grid  points  in  the  old  row. 

KN  current  number  of  grid  points  in  the  new  row. 


ITERO  maximum  number  of  iterations  for  the  integration  of  the  compatibility  relations 
(see  INVMAR  and  SEMINV;  also  used  in  RFUNC,  PLUMES). 

IM,  IP  search  indices  for  interpolation  subroutine  INTERP.  (see  INVMAR,  SEMINV). 

J  current  row  index  (also  index  of  a  marching  step). 

KF2  defined  as  2*  K  F;  not  used  in  present  version. 

IDEL,  JDEL  increments  for  printing  grid  point  I  and  row  J  (see  PRINT). 

JYXI  number  of  rows  to  be  printed  in  a  run. 

JXI  index  of  printing  row,  to  be  printed  next  (see  PRINT). 

ILEAD  index  I  at  the  first  grid  point  on  current  new  row,  where  the  SIMA  integration 

commences.  Initially  this  point  corresponds  to  the  leading  characteristic  of  the 
CRW.  (see  GRIDN,  BEGIN). 

ILEADF  value  of  ILEAD  for  current  old  row. 

KCLEAD  index  in  the  characteristic  array  for  the  characteristic  line  that  corresponds  to  the 
new  grid  point  I  *  ILEAD  (see  GRIDN).  Initially  KCLEAD  -  1. 


IM,  IP 


JYXI 


ILEAD 


ILEADF 

KCLEAD 


/STAG/  (INIDAT) 

RHOO,  NO  stagnation  density  and  number  density. 

PO,  TO,  AO  stagnation  pressure,  temperature  and  sound  speed. 

MDOT1  mass  ilow  rate  from  nng-nozzle  (only  from  the  x  >  0  half).  (See  PRINT). 
/ICHARA/  (BEGIN) 

KCHARP  number  of  C +  characteristic  lines  for  which  data  is  stored  (either  for  SIMA 
computation  or  for  subsequent  plotting). 

KCHARM  number  of  C  characteristic  lines  for  which  data  is  stored  (only  for  subsequent 
plotting). 

KCHARO  total  number  of  characteristics  for  which  data  is  stored,  i.e., 
KCHARO=  KCHARP+  KCHARM. 


mmmi 


3.4  Description  of  JET  subroutines 


MAIN  PROGRAM 

The  main  program  performs  two  functions.  The  first  section  (up  to  statement  1)  is  the  initial  set 
up;  it  is  performed  just  once.  The  second  section  is  the  marching  loop  with  the  step  index  J.  This 
program  can  be  read  as  a  flow  chart  of  the  overall  computational  procedure. 

INI  DAT  is  for  setting  up  run  data.  In  BEGIN  the  initial  conditions  for  the  marching 
computation  are  set  up.  A  single  marching  step  is  performed  by  calling  MARCH,  and  the  loading  of 
new  row  vectors  into  old  row  vectors  is  done  by  calling  MOVE.  The  call  to  YSTEP  is  for  the  first 
computed  marching  step.  All  remaining  calls  are  for  informative  tasks  (see  HMSET,  BREAK, 
OPACY,  PLUMES,  PRINT).  Run  is  terminated  when  either  YF.GE.YMAX  or  when  J.GE.JMAX. 

NOTE  ON  EXEC  :  The  only  special  feature  in  the  EXEC  is  retaining  the  output  unit  7  file  for 
optional  post-plotting.  The  printed  output  (unit  6)  is  the  system's  standard  (default). 


INIDAT 

Initial  data  definition  and  preliminary  data  computations.  The  data  is  defined  by  statements  rather 
than  by  reading  an  input  file.  The  meaning  of  major  parameters  was  described  in  Section  3.3  above. 
User  is  invited  to  modify  the  data  definitions,  particularly  of  run-control  parameters  such  as  YMAX. 
JYXI  and  YXI(JXI)  (for  printing  JYXI  selected  rows). 


BEGIN 

Here  all  initial  values  (prior  to  beginning  of  marching  schemes  integration)  are  loaded  into  all 
major  computational  arrays  (Section  3.1).  Also,  values  of  the  key  integer  parameters  KCHARP. 
KCHARM,  KCHARO,  ILEAD,  KCLEAD  and  KF  are  defined. 

In  the  first  loop  (loop  1)  we  define  an  initial  family  of  C+  characteristic  lines  for  the  lip-centered 
CRW,  bv  storing  the  Mach  number  of  the  Prandtl-Yleyer  fan  characteristics  in  the  array 


MCHARI(KC).  Note  that  the  fan  characteristics  are  generated  at  equal  RP  intervals,  since  the  flow 
variables  are  RM  and  RP.  However,  a  different  division  might  also  be  acceptable. 

The  next  step  is  the  definition  of  initial  values  for  all  characteristic  arrays,  first  the  C +  arrays 
(loop  2),  then  the  C  arrays  (loop  21).  The  C~  characteristic  lines  are  needed  just  for  informative 
output  (post-plotting),  so  the  present  version  contains  just  one  C  ~  line.  The  user  may  modify  that. 

The  remaining  grid  points  (altogether  KFO  grid  points  are  initially  available)  are  uniformly 
distributed  across  the  nozzle  opening,  and  the  row  arrays  are  loaded  with  the  corresponding  nozzle- 
exit  flow  variables  (loop  3). 

PRINT 

The  main  task  of  this  subroutine  is  the  printing  of  flow  variables  at  grid  points  of  selected  rows. 
The  printing  of  a  row  is  selected  when  YF  is  close  to  a  predefined  array  YXI(JXI).  Following  the 
printing,  JXI  is  updated  by  adding  1. 

For  comparison,  additional  flow  variables  are  printed  for  each  row.  These  are  computed  from  the 
analytic  approximation  to  a  ring-symmetric  CRW  (1J  ,  by  calling  MATCH.  Also,  lateral  molecular 
opacities  of  various  kinds  of  approximation  are  computed  by  calling  OPACX,  and  are  printed  for 
each  grid  point  within  the  CRW. 

Following  the  row  printing  (statement  120),  arrays  intended  for  post-processing  (plotting  of  special 
lines)  are  printed  and  subsequently  written  on  output  unit  7.  This  is  done  once  per  run.  just  before 
run  termination. 

FIN 

This  subroutine  is  called  when  an  error  is  encountered,  in  order  to  terminate  the  run.  Note  that 
the  run  is  terminated  by  deliberately  introducing  an  error  of  computing  SQRT(-l).  which  is  done  in 
order  to  trigger  the  printing  of  calling  sequences  by  the  operating  system. 


MARCH 


This  subroutine  performs  a  single  marching  step  by  calling  the  proper  computational  subroutines 
at  an  appropriate  sequence.  It  can  be  read  as  a  flow  chart  of  the  entire  computational  scheme.  First 
the  segment  of  the  new  row  suitable  for  SIMA  computation  is  calculated  by  calling  SEMIXV.  Then 
new  grid  points  for  that  part  of  the  new  row  for  which  inverse  marching  integration  is  to  be 
performed,  are  generated  by  calling  GRIDN.  The  results  of  the  SEMINV  computation,  which  were 
stored  in  characteristic  arrays,  are  now  loaded  into  row  arrays  by  calling  LOADC.  Finally,  the 
computation  of  the  new  row  is  completed  by  calling  INVMAR  which  computes  the  flow  at  the 
remaining  grid  points  by  the  inverse  marching  scheme. 


INVMAR 

This  is  one  of  the  two  central  subroutines  for  computing  the  flow  at  new  grid  points  (the  other  is 
SEMINV).  Here  the  inverse  marching  scheme  is  used.  The  computational  procedure  follows  the 
seven-step  description  given  in  Section  2.2  above.  Note  that  the  initial  value  of  the  search  indices  IM 
and  IP  is  not  redefmed  at  each  call  to  INTERP,  since  it  is  assumed  that  IM  and  IP  do  not  change 
much  at  consecutive  calls  to  INTERP,  so  that  search  efficiency  is  enhanced  by  not  starting  the  search 
from  an  arbitrary  point  (such  as  either  end  of  the  row). 


SEMINV 

This  is  the  subroutine  performing  the  SIMA  scheme  for  computing  the  flow  at  new  grid  points 
located  along  continuous  characteristic  lines  of  the  lip-centered  CRW  (at  prescribed  y-marching 
steps).  The  essence  of  the  computational  procedure  of  this  subroutine  was  given  as  a  seven-step 
description  in  Section  2.2.  The  same  remark  about  IM  given  in  the  preceding  INVMAR  description 
applies  here  as  well. 

The  main  loop  ( 100)  is  over  all  characteristic  lines,  including  some  C  lines  in  addition  to  the  C  + 
lines.  Thus,  the  array  CSIGNF(KC)  is  used  to  get  the  appropriate  expressions  for  either  C  +  or  C 
characteristics.  It  is  noted  that  while  normally  the  characteristic  segments  through  points  1  and  2  arc 
C ~  and  C  +  respectively,  this  is  reversed  when  aC~  rather  than  a  C  +  line  is  computed  via  the 
SIMA  scheme.  In  this  case,  which  is  characterized  by  having  CSIGNF(KC).LT.O,  the  Riemann 


invariants  integrated  along  segments  (i,4)  and  (2,4)  are  interchanged.  This  is  done  in  the  few 
statements  just  preceding  and  following  statement  21. 

An  additional  capability  of  this  subroutine  is  to  treat  a  change  of  a  C  +  characteristic  line  into  a 
C  line  upon  reflection  from  the  symmetry  plane  (x  *  0).  This  is  done  by  first  computing  a  new  grid 
point  having  X4.LT.0,  and  then  changing  its  sign  after  setting  CSIGN'N(KC)  =  —1  (statements  just 
preceding  statement  30).  It  is  also  possible  to  skip  the  computation  of  a  particular  characteristic  by 
setting  CSIGN'N(KC)  =  0.  This  feature  is  not  exploited  in  the  present  version. 

Finally,  we  note  that  not  all  characteristic  lines  computed  here  are  part  of  the  marching  flow 
computation.  Only  those  with  indices  KC  between  KCLEAD  and  KCHARP  are.  All  other 
characteristic  lines  are  computed  just  for  informative  purposes  (post-plotting). 

RFUNC 

Here  M,  MU,  TETA  are  computed  from  the  two  Riemann  invariants  RM,  RP.  The  computation 
of  M  is  performed  by  a  Newton-Raphson  iteration  using  Equations  (2.1-2)  and  (2.1-3)  given  in 
Section  2.1  above. 

INTERP 

This  subroutine  starts  by  finding  through  a  search  procedure  the  grid  interval  (I.  1  +  1)  that 
contains  a  given  point  X.  Then  the  Riemann  invariants  are  computed  for  this  point  by  linear 
interpolation,  and  returned  in  RM.  RP.  Note  that  X  may  be  negative,  which  accounts  for  the 
relatively  elaborate  search  logic  in  the  determination  of  I.  and  for  the  reflection  transformation  (as  in 
Eq.  (2.3-1)  above)  preceding  the  last  two  statements  of  the  subroutine. 

INTERX 

This  interpolation  routine  performs  an  inverse  task  to  that  of  INTERP.  in  that  it  finds  the  point 
X0  that  corresponds  to  a  given  linearly  interpolated  value  of  the  flow  variable  VARO.  It  is  used  in 
PLUMES  to  compute  the  location  of  a  breakdown  surface  point  on  a  new  row  of  x-centered  and  y- 
centered  grid  points 


BREAK 


i 

This  subroutine  computes  the  new  breakdown  parameter  array  BN(I).  The  computation  is  based 
on  the  description  given  in  Section  2.4  above. 
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i  OPACY 

;  Here  the  radial  (Y)  molecular  opacity  array  TH(J)  is  computed.  At  each  marching  step  J,  a  new 

boundary  grid  point  XTH(J)  is  added,  then  the  radial  opacities  at  all  preceding  boundary  points  are 
updated  by  adding  the  contribution  of  the  gas  layer  between  the  current  old  and  new  rows.  Note  that 
since  grid  points  on  adjacent  rows  are  not  located  on  equal-x  columns,  this  procedure  requires  x- 
!  interpolation  by  calling  INTERP. 

PLUMES 

This  is  a  user-defined  subroutine,  where  up  to  10  special  lines  can  be  computed  and  subsequently 
retained  on  output  unit  7  for  post-processing  (plotting).  The  type  of  the  line  ITYPL(IPL)  and  a 
parameter  VPL(IPL)  that  defines  the  line,  are  computed  through  user-inserted  statements  in  the 
section  preceding  statement  2000.  Then  an  additional  point  on  the  current  new  row  is  computed  for 
each  line  type.  The  available  types  are  clearly  stated  in  comments.  Note  that  characteristic  lines 
have  already  been  computed  in  SEMINV  using  the  SIMA  scheme,  regardless  of  whether  they  are  part 
of  the  solution  grid  to  the  flow  field,  or  are  just  computed  for  informative  purpose.  It  is  the  user  s 
choice  which  of  these  lines  (if  any)  are  to  be  saved  in  the  /PLUME/  arrays  for  subsequent  post¬ 
processing  (plotting). 


GRIDN 

This  subroutine  computes  the  grid  points  in  that  segment  of  the  new  row  for  which  the  flow  is 
computed  by  the  inverse  marching  scheme  (in  INVMAR).  Initially,  this  segment  extends  from  x  =  0 
to  the  new  row  grid  point  which  lies  on  the  leading  characteristic  of  the  lip-centered  CRW.  However, 
since  the  leading  characteristic  is  reflected  from  the  symmetry  plane  (x  =  0)  at  some  point,  this 
segment  steadily  shrinks  in  size  as  the  marching  proceeds.  The  remedy  is  to  declare  the  next-to-the- 


leading  characteristic  line  (KC~2)  as  the  beginning  of  the  segment  for  SIMA  integration,  by  setting 
KCLEAD  =  2.  This  process  of  increasing  KCLEAD  is  repeated  whenever  it  is  deemed  necessary. 
The  criterion  in  the  present  version  for  the  minimal  KCLEAD  is  that  the  inverse-marching  segment 
should  be  at  least  twice  DX1  -  the  average  CRW  grid  interval  (loop  1,  the  two  statements  following 
DX1  * ...).  Also,  I  LEAD  is  redefined  for  each  row  according  to  XLEAD/DX1  +  2  in  order  to  achieve 
a  row  of  relatively  uniform  grid  intervals  throughout.  The  result  is  that  the  number  of  grid  points  in 
a  row  is  initially  KFO,  but  eventually  it  decreases  due  to  both  increase  of  KCLEAD  and  decrease  of 
I  LEAD. 

YSTEP 

In  this  subroutine  the  next  marching  step  DYNEXT  is  computed  at  the  end  of  the  current 
marching  step.  It  is  defined  as  the  smallest  step  obtained  by  forward  intersection  of  C  ”  and  C  + 
characteristics  from  adjacent  grid  points.  Note  that  the  actual  value  of  DYNEXT  is  reduced  by  a 
"stability"  factor  STAB,  and  that  DY  is  also  limited  by  the  growth-rate  factor  DDY  and  by  DYMAX 
(see  MAIN  PROGRAM). 

MOVE 

Here  old  row  arrays  (loop  1)  and  old  characteristic  arrays  (loop  2)  are  loaded  with  values  of  flow 
variables  from  corresponding  new  arrays,  in  preparation  for  the  next  marching  step.  As  a  result  of 
this  organizational  feature,  informative  computations  (e.g.  BREAK.  OPACY)  that  require  both  new 
and  old  rows,  have  to  be  performed  prior  to  calling  MOVE. 

OPACX 

Here  lateral  (X)  opacities  that  correspond  to  the  number  of  expected  collisions  of  a  fast  molecule 
invading  the  CRW  in  the  -X  direction,  are  computed.  All  opacities,  except  XIAPP(I).  are  computed 
by  numerical  integration.  In  loop  1  we  compute  the  opacity  contribution  of  the  segment  lying  just 
outside  the  computational  boundary  characteristic  (MFIN),  assuming  a  Prandtl-Meyer  flow.  This 
additional  opacity  is  denoted  XIO.  If  the  flow'  is  ring-symmetric,  XIO  is  recalculated  using  the  analytic 
approximation  (1)  to  estimate  the  flow  field  at  the  fringes  of  the  ring-symmetric  CRW  (see  also  the 
closed  form  expression  for  t  in  [1]  ). 


The  computation  of  opacity  arrays  starts  after  statement  14.  First,  the  opacity  at  each  grid  point 
is  set  to  XIO.  Thus,  even  though  the  numerical  flow  computation  does  not  include  the  fluid  outside 
the  boundary  characteristic  line,  the  opacity  integration  includes  an  estimate  of  that  “missing"  part, 
i.e.,  of  XIO.  In  typical  case  computations  of  a  ring- symmetric  CRW  [1J  we  found  that  the  maximum 
value  of  XIO  was  about  0.16.,  w’hich  indicated  that  as  far  as  interaction  with  invading  ambient 
molecules  is  concerned,  the  approximation  MFIN’  =  34  was  a  reasonable  substitute  for  MFIN=  oo. 

The  next  step  is  the  computation  by  numerical  integration  of  three  approximations  to  the  lateral 
opacity  :  XI(I,JXI),  XIPM(ITXI),  XIGRP(I,JXI).  (Note  that  when  the  flow  is  ring-symmetric,  the 
approximation  XIPM(IrJXI)  obtained  by  assuming  a  Prandtl-Mever  flow  is  usually  grossly 
exaggerated).  The  opacity  XIGRP(I,JXI)  is  based  on  the  analytic  approximation  to  a  ring-symmetric 
CRW  (1]  ,  and  is  reasonably  close  to  XI(I,JXI)  which  is  obtained  from  the  numerical  solution  to  the 
flow  field.  Finally,  a  simplified  closed-form  integration  of  lateral  opacity  [1]  is  computed  as 
XIAPP(I,JXI)  (loop  3).  Thus,  the  quantitative  difference  between  XI(I,JXI)  and  XIGRP(I,JXI)  is  an 
indication  to  the  degree  of  accuracy  achieved  by  the  analytic  approximation  to  a  ring-symmetric 
CRW  [1]  ,  while  the  difference  between  XIGRP(I,JXI)  and  XIAPP(IrIXI)  indicates  the  level  of  error 
introduced  by  the  closed-form  integration  of  lateral  opacity  [1]  . 

LOADC 

Here  flow  variables  of  new  grid  points  computed  via  the  SIMA  scheme  (SEMINV)  are  loaded  into 
new  row  arrays  from  corresponding  characteristic  arrays. 

NUFUNC 

This  function  computes  the  modified  v(M)  value  as  given  by  Eq.  (2.1-2).  Note  that  presently 
NL'0-0  (see  INI  DAT). 

HMSET 

This  subroutine  is  called  just  once  from  the  MAIN  PROGRAM.  Its  task  is  to  set  up  the  arrays  in 
/ GRP/,  so  that  the  function  H(M)  (1J  can  be  evaluated  by  interpolation  (in  HINTER).  There  is  also 
an  informative  printout  of  various  derivatives  (see  Ch.  3  of  (1)  )  generated  in  this  subroutine. 


MFUNC 


This  subroutine  is  called  by  HMSET  in  order  to  compute  functions  of  Mach  number  that  serve  in 
the  computation  of  H(M).  The  output  variable  F  is  the  integrand  for  the  integration  leading  to 
H(M). 

HINTER 

This  subroutine  computes  H(M)  by  linear  interpolation  in  inverse  Mach  number,  using  the  /GRP/ 
arrays  computed  in  HMSET. 


MATCH 

This  subroutine  is  called  from  PRINT  to  compute  the  Mach  number  according  to  the  analytic 
approximation  of  a  ring- symmetric  CRW  flj  ,  for  point  (YF,XF(I)).  MOB  is  the  associate  Mach 
number  M(0,p),  which  is  preserved  in  the  array  MCHARI(KC)  for  all  CRW  characteristics  that  are 
used  in  the  SIMA  computation.  Hence  the  Mach  number  M(a,P),  denoted  by  MAB  can  be 
computed  directly  from  the  analytic  approximation  [1]  to  the  area  function  at  (YF.XF(I))  by  calling 
AREAF.  Since  typically  M(O.P)  is  not  known,  we  also  compute  the  Mach  number  via  the  inverse- 
problem  procedure  Ill  ,  denoting  the  resulting  Mach  numbers  by  suffix  I  :  M0B1  for  M(0,p)  and 
MAB  I  for  M(a,p).  The  inverse-problem  iterative  procedure  [1]  is  performed  in  loop  1,  resulting  in 
MOBI.  From  MOBI  the  value  of  MABI  is  computed  through  the  area  function  approximation  as  for 
MAB  above. 


AREAF 

This  subroutine  computes  the  Mach  number  M  that  corresponds  to  the  area  function  F  (Eq. 
(3.2-1)  of  {1]  ).  The  computation  is  done  by  Newton-Raphson  iterations,  and  it  has  been  found  to 
converge  when  M.GT.l  (and  when  M  -  1  is  not  much  smaller  than  1). 


4.  THE  JET  CODE  LISTING 

CBOPTIONS  LIST 
C  J CTO 18 

S  SEMI-INVERSE  MARCHING  CHARACTERISTICS  METHOD  FOR  RING  JETS 

C  VARIABLe|MANN  INVARIANTS  RM*(NU+TETA),  RP.(NU-TETA)  AS  FIELD 

IMPLICIT  REAL*8<A-H,L-Z,$) 

REALK4  XPL,YPL 

COMMON  /PLUME/XPL ( 1002, 10),  YPK1002) 

COMMON  /IPLUME/KPL,ITYPL(10) 

COMMON  /VECS/XF< 101 ) , RMF( 101 ) , RPFC 101 ) ,MFC 101 ) ,MUF( 101 ) , 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3  TETAN(101),BNC101),XTEMP(101) 

COMMON/THICKY/XTH( 1002), TH< 1002) 

REAL*?,  yxi,xi,xipm,xigrp,xiapp,xif 

COMMON  /THICKX/YXK20) ,XI(101,20),XIPM(101, 20) , XIGRP(  101 ,20) 

1  ,  XIAPPdOl,  20),XIF(101,20) 

COMMON  /GAMA/G, G1 ,  G2, G3, G4, G5, G6 ,G7,G8,G9,G10,G11,G12,G13,G14,G15 

1  G16 , G17 ,G1$ ,G19 ,G2Q 

COMMON  /PAR/PAI,PAI2, DEG,XC,YC,MEXIT,MFIN, YMAX, DYO, DY, DYNEXT, 

2  ?ETSm  TEUl5;DD??DTn«TA1  • SI0H4'  FR‘C0- EPSU '  ^ 
COMMON  /STAG/RHOO , NO , PO , TO , AO , MDOT1 

COMMON  /IPAR/ JMAX, KFO,ITERO,KF,KN,IM,IP,J, 

^COMMON  ,R<lH,YFFYK!ExF:BXNL'JYXI'JXI'UEAD'IlEADF',CCtEAI> 

COMMON  /CHARAC/XCHARFC92) , YCHARF( 92) , XCHARNC 92) , YCHARNC 92) , 

\  RMCARF(92),RPCARFC92),RMCARN(92),RPCARN(92), 

I  '  TCHARN(  92 ) ,  MUCARF(  92)  , MUCARN(  92 ) , 

l  MCHARI ( 9| J ' CSIONF  i C  92 ) , MCHARN( 9  2 ) , MCHARFC  92 ) , 

COMMON  /ICHARA/KCHARP , KCHARM, KCHARO 
COMMON  /GRP/DMINV , MHINV(lOl) , HMV( 101 ) 

COMMON  /IGRP/KHM 
C 

PRINT  101 
101  FORMAT ( 1 1  * ) 

J  =  1 

IF(J.EQ.l)  STOP 
CALL  INIDAT 
PRINT  101 
CALL  HMSET 
PRINT  IOI 
CALL  BEGIN 
CALL  MARCH 
CALL  OPACY 
CALL  PLUMES 
CALL  PRINT 
J=2 

CALL  PLUMES 
CALL  MOVE 
CALL  OPACY 
CALL  PRINT 
CALL  YSTEP 
1  J=J+1 

C  DY  WAS  DETERMINED  BY  THE  PREVIOUS  CALL  TO  GRIDN. 

„  t.,t5YbDMIN1  (  DYNEXT ,  DY*DDY ,  DYMAX) 

C  INTEGRATE  BY  ONE  Y-STEP 
CALL  MARCH 

C  BREAKDOWN  PARAMETER  ( BF( I) ) . 

CALL  BREAK 

C  SPECIALLY  DESIGNATED  LINES  (FOR  PLOTTING). 

CALL  PLUMES 

C  STORE  NEW  LINE  (N)  IN  OLD  LINE  (F). 

CALL  MOVE 

C  COMPUTE  RADIAL  MOLECULAR  OPACITIES 
CALL  OPACY 

C  Y_SIF(YFSGEAYMAX)EJMAX=jMAX  IS  USED  AS  END"0F~RUN  CRITERION. 

C  PRINT  FIELD  AT  MOST  RECENT  Y. 

CALL  PRINT 
C  NEXT  Y-STEP. 


JET0001 
JET0002 
.  JET00Q3 
JET0004 
JET0005 
JET0006 
JET0007 
JET0008 
JET0009 
JET0010 
JETQOll 
JET0012 
JET 001 3 
JET0014 
JET0015 
JET0016 
JET0017 
i,  JET0018 
JET0019 
JET0020 
JET0021 
JET0022 
JET0023 
JET0024 
JET0025 
JET0026 
JET0027 
JET0028 
JET0029 
JET0030 
JET0031 
JET0032 
JET0033 
JET0034 
JET0035 
JET 00 36 
JET 0037 
JET0038 
JET0039 
JET0040 
JET0041 
JET0042 
JET0043 
JET0044 
JET0045 
JET0046 
JET0047 
JET0048 
JET0049 
JET0050 
JET0051 
JET0052 
JET0053 
JET0054 
JET0055 
JET0056 
JET0057 
JET0058 
JET0059 
JET0060 
JET0061 
JET0062 
JET0063 
JET0064 
JET0065 
JET0066 
JET0067 
JET0068 
JET0069 
JET0070 
JET0071 
JET0072 


FILE.  JETPR 


FORTRAN  A1 


hv* 


A-V 

It? 


KY 


C 

c 


CALL  YSTEP 

IF(J.LT.JMAX)  GO  TO  1 
STOP 


-END 


iMJWT' 


JET0073 

JET0074 

JET0075 


SUBROUTINE  INIDAT 
SUBROUTINE  NUMBER  1 

IMPLICIT  REAL*8CA-H,L-Z,$) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(1Q1),MUFC101), 

1  TETAFC101),BF(101), 

2  XN( 101 ) , RMNC 101 ) , RPN( 1 01 ) ,MN( 101 ) ,MUN( 101 ) , 

3  TETAN(101),BN(101), XTEMPC 101 ) 

COMMON/THICKY/XTH( 1002) ,TH(1002) 

REAL*4  YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON  /THICKX/YXI ( 20 ) , XI ( 101 , 20 ) , XIPM( 101 , 20 ) , XIGRPCl 01 , 20 ) 

1  ,XIAPP(101,20),XIF(101,20) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2, DEG, XC, YC, MEXIT, MFIN, YMAX, DYO, DY, DYNEXT, 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG, EPSIL,NUO, 

2  TETSYM, TETL IM, DDY, DYMAX 
COMMON  /STAG/RHOO, NO, PO, TO, AO , MD0T1 
COMMON  /IPAR/JMAX, KFO, ITERO , KF, KN, IM, IP, J , 

1  KF2, 1  DEL , JDEL , JYXI , JXI , I  LEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF, DXN 


JJFJAtnf, 


PAI=4 . DOXDATANC 1 .DO) 

PAI2=2.D0*DATAN(1 .DO) 

DEG=180 . DO/PAI 
AR=8 . 3143D3 
AV=6 . 022D  26 
AW=7 . 27D0 
RHO0=0 . 0075D0 
T0=2300 . DO 
G=1 . 54D0 
D=2 . 5D-10 
MEXIT=4 . DO 
MFIN=34 . DO 
XC=0 . 5D0 
YC=2 . 5D0 

DELTA=0  CORRESPONDS  TO  PLANE  SYMMETRY 

DELTA=1  CORRESPONDS  TO  CYLINDRICAL  SYMMETRY 
DELTA=1 . DO 


FRACG=0 . 6D0 
EPSIL=l.D-8 
ITER0=20 
KF0=1 01 
JMAX=1001 
STAB=0 . 50D0 
DDY=1 . 05D0 
DYMAX =0 . 5D0 
YMAX=50 . DO 
DYO=YC/250.DO 
I  DEL  =  1 
JDEL=1 

POINTS  FOR  PRINTING  FLOW  FIELD  AT  YF=YXI(JXI) 

JXI  =  1 
JYXI  = 1 1 
DYXI =5 . DO 
YXI ( 1 ) =YC+0 . 5D0 
YXI ( 2 ) =YXI C 1 ) +2 . DO 
10  =  2 

DO  1  1=10, JYXI 

YXI(I)=YXI(IO ) +DYXI XDFLOAT  (  I-1 0 ) 

CONTINUE 

IFCKFO.GT.lOl)  CALL  FIN(lOl) 

I F( JMAX . GT . 1 00 1 )  CALL  FINC102) 

I F( FRACG . GT . 1 . DO  .OR.  FRACG. LT.O.)  CALL  FINC103) 
IF< JYXI .GT.20)  CALL  FINC104) 

IF(DELTA*<1. DO-DELTA). HE. 0. )  CALL  FIN(105) 

NO=RHOO#AV/AW 

AO=DSQRT(G*AR*TO/AW) 

PO=AR*RHOO*TO/AW 


21 


JET0077 
JET0078 
JET0079 
JET0080 
JET 0081 
JET0082 
JET0083 
JET0084 
JET0085 
JET0086 
JET0087 
,  G15, JET0088 
JET0089 
JET0090 
JET0091 
JET0092 
JET0093 
JET0094 
JET0095 
JET0096 
JET0097 
JET0098 
JET0099 
JET0100 
JET0101 
JET0102 
JET0103 
JET0104 
JET0105 
JET0106 
JET0107 
JET0108 
JET0109 
JET0110 
JET0111 
JET0112 
JET0113 
JET0114 
JET0115 
JET0116 
JET0117 
JET0118 
JET0119 
JET0120 
JET0121 
JET0122 
JET0I23 
JET0124 
JET0125 
JET0126 
JET0127 
JET0128 
JET0129 
JET0130 
JET0131 
JET0132 
JET0133 
JET0134 
JET0135 
JET0136 
JET0137 
JET0138 
JET0139 
JET0140 
JET0141 
JET0142 
JET0143 
JET0144 


.'■■.v.  v.  /.  •* 


JETPR 


FORTRAN  A1 


C 

21 

22 


23 


SIGMA*PAIXDXX2 

LAMDAO*l.DO/(DSQRT(2.DO)XSI6MAXHO) 

Gl*(0-1 . 00 )/2 . DO 
G2«(G+1,D0)/(2.D0X(G-1.D0)) 

63*0/2 . DO 

G4*(G+1 . D0)/(G-1 . DO) 

G5*DSQRT( (G+l . D0)/(G-1 . DO) ) 

G6*l . D0/(G“1 . DO) 

G7*2 . DO/ (G+l . DO) 

G8=(0.5D0x(G+l.D0)xx2/CG-l.D0))xx(i.D0/CG+l.D0))x 
1  ( (G+l . D0)/(G-1 . DO) )xx( CG-1 . DO )/ (G+l . DO ) ) 

G9*(G+3 . DO )/(2 . D0X(G-1 .DO)) 

G10*(7 . DO-3 . DOXG)/( 2 . D0X(G-1 . DO)  ) 

Gil =( 2 . DO/ (G+l . DO))xx(i . D0/(G-1 .DO)) 

G12*DSQRT ( (G+l . D0)/(G-1 . DO) )-l . DO 
G13=(2 . D0-G)/(2 . DOx(G-l . DO) ) 

G14=G/(2 . D0X(G-1 . DO) ) 

G15=(G+1 . D0)/(3 . DO-G) 

G16*(G+1 . D0)/4 . DO 
G20=LAMDA0XDSQRT( PAIXG/8 .DO) 

ZETA1=G5*DATAN(  DSQRKMEXITXX2-1 .  DO  )/G5) 
AMU1=DARSIN( 1 . DO/MEXIT) 

PSI1=PAI2+AMU1 

ZETAF=G5XDATAN( DSQRT(MFINXX2-1 . DO )/G5) 

PSIF=PSI1+ZETA1-ZETAF 

NUO-O. 

TETLIM=NUFUNC(MEXIT)+PAI2-NU0 
PSILIM*TETLIM 
TETSYM=PAI-2 . DOXTETLIM 
G0REM=1 . D0+G1XMEXITXX2 
RH01 -RHOO/GOREMKXG6 
V1=MEXITXAO/DSQRT(GOREM) 

P1=PO/GOREMXX(G/(G-1 . DO ) ) 

T1=TO/DSQRT(GOREM) 

YYC=2 . DOXPAIXYC 
IF(DELTA . EQ . 0 . )  YYC=1.D0 
MD0T1=YYC*XC*RH01*V1 


JET0145 
JET0146 
JET 01 47 
JET 01 48 
JET0149 
JET0150 
JET0151 
JET 01 52 
JET0153 
JET0154 
JET0155 
JET0156 
JET0157 
JET 01 58 
JET0159 
JET0160 
JET0161 
JET0162 
JET0163 
JET0164 
JET0165 
JET0166 
JET0167 
JET0168 
JET0169 
JET0170 
JET0171 
JET0172 
JET0173 
JET0174 
JET0175 
JET0176 
JET0177 
JET0178 
JET0179 
JET0180 
JET0181 


JET0182 


PRINT  21,AR,AV,AW,G,RHOO,NO,PO,TO,AO,D  JET0183 

F0RMAT(/1X, 'THERMODYNAMIC  DATA.'/  JET0184 

1  IX, 'AR,AV,AH,G=',2X,2D14.5,2F9.3/  JET0185 

2  IX, 'RHOO,NO,PO,TO,AO,D=',6D13.5)  JET0186 

PRINT  22,XC,YC,MEXIT,RH01,P1,T1,V1,MD0T1,PSI1XDEG,PSIFXDEG,  JET0187 

1  PSILIMXDEG  JET0188 

F0RMAT(/1X, 'CORNER  DATA.  XC, YC= ' , 2F9 . 2/  JET0189 

1  IX, 'EXIT  CONDITIONS: ',  JET0190 

2  2X, 'MEXIT,RH01, PI, Tl, VI, MD0T1=',F9. 3,5013.4/  JET0191 

3  IX, 'CENTERED  FAN  LIMITS:',  JET0192 

4  2X, 'PSI1,PSIF,PSILIM=',3F10.3)  JET0193 

PRINT  23, DELTA, KFO , JMAX, ITERO , DYO, YMAX, STAB, DDY  JET0194 

F0RMAT(/1X, 'INTEGRATION  DATA.  SYMMETRY  INDEX:  DELTA- ', F4 . 1/  JET0195 

1  IX, 'NUMBER  OF  POINTS  IN  X  AND  Y  DIRECTIONS.  KFO,JMAX=',  JET0196 

2  215/  JET0197 


RETURN 

END 


IX, 'MAX.  NUM.  OF  ITERATIONS  ITER0=»,I5/ 

IX, 'INITIAL  Y-STEP  AND  MAXIMUM  Y:  DYO , YMAX= • , 2D14 . 5/ 
IX, 'Y-STEP  STABILITY  FACTORS  STAB, DDY= • , 2F7 . 3) 


SUBROUTINE  NUMBER  2 

IMPLICIT  REAL*8(A-H,L-Z,$) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1  TETAF(101),BF(101), 

2  XN(101  ),RMN(101 ),RPN(101 ) , MN( 101 ) , MUN( 101 ) , 

3  TETAN(101),BN(101), XTEMP( 101) 

COMMON/THICKY/XTH( 1002) ,TH( 1002) 

COMMON  /GAMA/G, G1 ,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PA I , PAI 2 , DEG, XC, YC, MEXIT,MFIN, YMAX, DYO , DY, DYNEXT, 

1  STAB, DELTA, PSI1,PSIF»ZETA1, SIGMA, FRACG, EPSIL,NUO, 

2  TETSYM, TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO, NO,PO,TO,AO,MDOT1 
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COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,  J,  JET0217 

1  KF2, I DEL, JOEL , JYXI , JXI ,  XLEAD, ILEADF, KCLEAD  JET0218 

COMMON  /ROW/YF,YN, DXF, DXN  JET0219 

COMMON  /CHARAC/XCHARFC92),YCHARFC92),XCHARNC92),YCHARNC92),  JET0220 

1  RMCARF( 92) / RPCARFC  92) ,  RMCARNC  92) , RPCARN( 92 ) ,  JET0221 

2  TCHARFC92) » TCHARN(92) ,MUCARFC92) ,MUCARNC92) ,  JET0222 

3  CSIGNNC  92) , CSIGNFC  92) , MCHARNC92) , MCHARFC92) ,  JET0223 

4  MCHARK92)  JET0224 

COMMON  /ICHARA/KCHARP, KCHARM, KCHARO  JET0225 

C  JET0226 

C  DEFINE  INITIAL  CHARACTERISTIC  PARAMETERS.  USE  INTERPOLATION  OF  JET0227 

C  RIEMANN  INVARIANT  ACROSS  THE  FAN.  JET0228 

KCHARP* I DINT ( FRACG*DFLOATCKFO_l)+l . D-6)+l  JET0229 

KCHARO =KCHARP+1  JET0230 

KCHARM-KCHARO-KCHARP  JET0231 

IF< KCHARP.LT. 2  )  CALL  FINC200)  JET0232 

I F( KCHARO . GT . 92)  CALL  FINC210)  JET0233 

IFCKCHARM.LT.  1)  CALL  FINC205)  JET0234 

NU1*NUFUNCCMEXIT)  JET0235 

RMl=NUO  JET0236 

TET1=RM1-NU1  JET0237 

RP0=NU1-TET1  JET0238 

NUFIN=NUFUNCCMFIN)  JET0239 

RPFIN=NUFIN-CRM1-NUFIN)  JET0240 

DRP=C  RPFIN-RPO )/DFLOAT ( KCHARP-1 )  JET0241 

DO  1  KC=1, KCHARP  JET0242 

RPl=RPO+DRP*DFLOATCKC-l)  JET0243 

CALL  RFUNCC RM1 , RP1 , Ml < MU1 *  TETA1 )  JET0244 

MCHARI CKC)=M1  JET0245 

1  CONTINUE  JET0246 

C  DATA  FOR  C+  CHARACTERISTICS.  JET0247 

C  THE  RIEMANN  INVARIANTS  ARE  DEFINED  IN  SUCH  A  WAY  THAT  BOTH  VANISH  AT  JET0248 
C  INFINITE  MACH  NUMBER.  JET0249 

RM1-NUO  JET02S0 

DO  2  KC=1 , KCHARP  JET0251 

CSIGNFC KC)=1 . DO  JET0252 

XCHARFCKC)=XC  JET0253 

YCHARF(KC)=YC  JET0254 

IFCMCHARICKC) . EQ . 0 . )  CALL  FINC231)  JET0255 

NU=NUFUNC( MCHARI ( KC) )  JET0256 

TET =RM1-NU  JET0257 

RP1=NU-TET  JET0258 

CALL  RFUNCC RM1 , RP1 , Ml , MU1 » TETA1 )  JET0259 

MCHARFCKC)=M1  JET0260 

MUCARFCKC)=MU1  JET0261 

TCHARFCKC)=TETA1  JET0262 

RMCARF(KC)=RM1  JET0263 

RPCARFCKC)=RP1  JET0264 

2  CONTINUE  JET0265 

C  DATA  FOR  C-  CHARACTERISTICS.  JET0266 

KC1=KCHARP+1  JET0267 

XCHARFCKC1 )=0 . 8D0XXC  JET0268 

DO  21  KC=KC1, KCHARO  JET0269 

CSIGNFCKC)=-1 .DO  JET0270 

MCHARI CKC) =MEXIT  JET0271 

MUCARF(KC)=DARSINC1 .  DO/MCHARICKC) )  JET0272 

TCHARFCKC)=PAI2  JET0273 

YCHARFCKC)=YC  JET0274 

MCHARFCKC) =MEXIT  JET0275 

RMCARFCKC)=RM1  JET0276 

RPCARFC KC)=NUFUNCCMEXIT)-CTCHARFCKC)-TETLIM)  JET0277 

21  CONTINUE  JET0278 

C  DEFINE  GRID  AND  INITIAL  CONDITIONS  AT  EXIT  PLANE.  JET0279 

KFAN=KCHARP-1  JET0280 

ILEAD=KFO-KFAN  JET0281 

KCIEAD=1  JET0282 

KF=KF0  JET0283 

KF2=2*KF  JET0284 

YF=YC  JET0285 

DO  3  1=1, KF  JET0286 

KC=KCLEAD+I-ILEAD  JET0287 

IFCKC.GT. KCHARP)  CALL  FINC241)  JET0288 


23 
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IF(KC.GE.l)  60  TO  31 

XF<I)*DFL0ATCI-1)*XC/DFL0ATCILEAD-1) 

MF(I)*MEXIT 

TETAF<I)=PAI2 

GO  TO  32 

CONTINUE 

XF(I)=XC 

MFCI)=MCHARF(KC) 

TETAFCI)=TCHARFCKC) 

CONTINUE 

RMF( I )aNUFUNC<MF( I ) )+(TETAF< I )-TETLIM) 
RPF( I )aNUFUNC(MF( I ) )-(TETAF( I )-TETLIM) 
MUF( I ) *DARSIN( 1 . DO/MF< I ) ) 

BF(I)aO. 

CONTINUE 

DYaDYO 

DO  4  KCal ,  KCHARO 
CSIGNN(KC)aCSIGNF(KC) 

CONTINUE 
DO  5  1*1, KN 
BN(I)=0 . 

CONTINUE 

RETURN 

END 


PRINT- 


SUBROUTINE  NUMBER  3 

IMPLICIT  REAL*8(A-H,L-Z,$) 

REALX4  XPL,YPL 

COMMON  /PLUME/XPK1002, 10),YPL(1002) 

COMMON  /IPLUME/KPL, ITYPL( 10) 

COMMON  /VECS/XF(101) , RMF( 101 ) ,RPF(101),MF(101),MUF(101), 

1  TETAF(101),BF(101), 

2  XN( 101) , RMN( 101 ),RPN(101), MN( 101) ,MUN( 101) , 

3  TETAN(1Q1),BN(101),XTEMP<101) 

COMMON/THICXY/XTH( 1002), THC 1002) 

REAL *4  YXI , XI , XIPM, XIGRP , XIAPP , XI F 

COMMON  /THICKX/YXI(20),XI(101,20),X1PMC101,20),XIGRP(101,20) 

1  ,XIAPP(101,20),XIFC101,20) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14 
1  G16 , G17 , G18 , G19 , G20 

COMMON  /PAR/PA I , PAI2, DEG, XC, YC, MEXIT » MFIN, YMAX, DYO, DY, DYNEXT , 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO , 

2  TETSYM, TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO , NO , PO , TO , AO , MD0T1 

COMMON  /CHARAC/XCHARFC92) , YCHARFC92) , XCHARN( 92 ) , YCHARNC 92) , 

1  RMCARF( 92 ) , RPCARF( 92 ) , RMCARNC  92) , RPCARNC  92 ) , 

2  TCHARF( 92),TCHARN(92), MUCARFC  92) , MUCARNC  92 ) , 

3  CSIGNN( 92 ),CSIGNFC92), MCHARNt  92 ) , MCHARF( 92 ) , 

4  MCHARK92) 

COMMON  /IPAR/JMAX, KFO , ITERO , XF , KN , IM, IP, J, 

1  KF2, 1  DEL , JDEL ,JYXI,JXI,ILEAD, IL  EADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF , DXN 

SUM=0. 

KF1 =KF-1 

DO  10  1=1, KF1 

DX=XF(I+1)-XF(I) 

G0REM=1 .D0+G1*MFCI)X*2 
G0REM1=1 .D0+G1XMF(I+1)X*2 

RATEM=RHOOXAOXMFCI  ) XDSIN( TETAF( I  ) )/GOREM  *X(G6  +  0.5D0) 
RATEP=RHOO*AOXMF(I+1)XDSIN(TETAFCI+1))/GOREM1XX(G6+0.5DO) 
SUM=SUM+DXX(RATEM+RATEP)/2 . DO 
)  CONTINUE 

YYF=2.D0XPAI*YF 
IFCDELTA . EQ . 0 . )  YYF=1.D0 
MD0TFR=YYFXSUM/MD0T1 

PRINT  11,  J, KCLEAD, KF, ILEAD, YF, DY, XF( KF) , MFC KF) , MDOTFR 
L  FORMAT (1X,'J,KCLEAD,KF,ILEAD,YF,DY,XF(KF), MBOUND,MDOTR= ' , 

1  415, 5D12 . 4 ) 

PRINT  FLOW  FIELD  AT  Y=YF 
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IFCJ.EQ.JMAX)  JXI=MIN0C JXI , JYXI ) 

IFCJ.EQ.l  .OR.  J.EQ.JMAX)  GO  TO  121 
IF( JXI. GT. JYXI)  GO  TO  120 
IF(YXKJXI)  .GT.YF+0.5D0XDY)  GO  TO  120 
121  CONTINUE 

YXI( JXI )SYF 
CALL  OPACX 

C  COMPUTE  MACH  NUMBER  FOR  CYLINDRICAL  EXPANSION  MCYL . 

F= C YF/YC ) * ( G7  * ( 1 . D0+G1*MEXIT**2) )x*G2/MEXIT 
CALL  AREAFCF,MCYL) 

PRINT  22, JXI ,  KCLEAD, ILEAD, KF, MCYL , YF 
22  FORMATC/1X, 'PRINTING  NUMBER  JXI , KCLEAD, ILEAD, KF=‘ , 414, 
1  5X, 'MCYL , YF= ' , 2D14 . 5/) 

PRINT  1 


FORMAT (/IX, • 

I  ','  KC  ',' 

XFC I ) 

TETAFCI) 

1 

9 

1 

« 

MFCI) 

MAB 

*  9 

2 

f 

MABI 

MOBI 

1 

; 

3 

V 

XICI) 

XIGRPC I ) 

*  9 

4 

V 

XIAPPCI) 

XIPMCI) 

•/) 

IDEL1=IDEL 

IFCJ.EQ.l. OR. J.EQ.JMAX)  IDEL1*1 
DO  20  1=1 , KF, IDEL1 
KC=KCLEAD+C I-ILEAD) 

IFLKC.LT. KCLEAD)  KC=0 

M0B=1 . DIO 

M0BI=1 . DIO 

MAB=1 . DIO 

MABI=1 . DIO 

MPM=MF< I ) 

IFCKC.EQ.O)  GO  TO  23 
MOB=MCHARI(KC) 

IF(J.EQ.l)  GO  TO  23 
PSIPM=PAI2-DATAN( (XF( I )-XC)/(YF-YC) ) 

ZETA=PSI1+2ETA1-PSIPM 
MPM=DSQRTCCG5*DTANCZETA/G5))*X2+1 .DO) 

CALL  MATCH(I,MOB,MAB,MOBI,MABI ) 

23  CONTINUE 

PRINT  21,I,KC,XF(I), TETAFC I )#DEG»  MFC I ), MAB, MABI , MOBI , 

1  XICI.JXI) , XIGRPC I , JXI ) , XIAPPC I,JXI),XIPMCI,JXI) 

21  FORMATC1X, 214,10012. 4) 

20  CONTINUE 

IFCJ.EQ.l)  GO  TO  120 
IF(J.EQ.JMAX)  GO  TO  120 
JXI=JXI+1 
120  CONTINUE 

IF(J.LT.JMAX)  GO  TO  200 
PRINT  101 

101  FORMAT ( ' 1 '  ) 

PRINT  102 

102  F0RMATC1X, 'RADIAL  MOLECULAR  THICKNESS  J , XTHC  J  ) ,  THC  J  )  =  '/ ) 
PRINT  202, (J J, XTHC JJ ), THC J J  )  , J J  =  1 , JMAX ) 

202  F0RMATC/5CI5, Dll. 4,010.3)) 

PRINT  101 

PRINT  103, ( I PL , ITYPL (IPL),IPL=1,KPL) 

103  FORMATC IX, 'PLUME  TYPES  IPL , ITYPL C IPL )=' , 

1  2C/1X, 5C  5X, 214 )) ) 

PRINT  104 

104  FORMATCIX, 'PLUME  POINTS  J , YPL C J ), XPL C J , 1 ), XPL ( J , 2 ),...  =  •/ ) 
JDEL1=1 

DO  203  JJ=1, JMAX, JDEL1 

PRINT  204, JJ,YPL(JJ), C XPL C JJ, IPL ) , IPL  =  1 , KPL  ) 

204  FORMAT (1X,I5,2X,E12.4,10E11.3) 

203  CONTINUE 

C  WRITE  ON  TAPE7  FOR  SUBSEQUENT  PLOTTING. 

C  NO  MORE  THAN  80  CHARACTERS  PER  LINE  ON  TAPE7  . 

WRITEC 7 , 205 )  JMAX, KPL 

205  FORMATC 81 10/81 10) 

WRITEC 7, 205)  C ITYPL C IPL ), IPL =1 , KPL  ) 

DO  210  JJ=1 , JMAX 

WRITEC 7, 211)  YPLCJJ),CXPLCJJ,IPL),IPL=1,KPL) 

211  FORMAT  C6E13.6/2X,6E13.6/2X,6E13.6/2X,6E13.6) 

210  CONTINUE 
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226 


JXI0,KF0=',2I8) 


221 

227 


WRITE  LATERAL  <X)  OPACITIES 
JXIO=JXI 

WRITEC7 , 205)  JXIO,KFO 
PRINT  226 ,  JXIO,KFO 
F0RMAT(///1X, 'LATERAL  (X)  OPACITIES 
DO  220  JXI=1,JXI0 
WRITE(7, 221 )  JXI,YXI(JXI) 

FORMAT (110, El 3. 6) 

PRINT  227,  JXI.YXIUXI) 

F0RMATC//1X, • JXI , YXI ( JXI ) = ' , 18 , E15 . 6/ ) 

DO  225  1=1, KFO 

WRITEC7, 211 )  XI F( I, JXI) , XI(I, JXI),XIPM(I, JXI),XIGRP(I, JXI), 
XIAPP( I , JXI ) 


225 

220 

200 


PRINT  211,  XIF(I , JXI ) , XI ( I, JXI ) , XIPM( I , JXI ) ,XIGRP( I , JXI ) , 
1  XIAPP( I , JXI ) 

CONTINUE 
CONTINUE 
CONTINUE 
RETURN 
END 
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SUBROUTINE  FTNiTFIN) - 

SUBROUTINE  NUMBER  4 

STOP  WHEN  ERROR  IS  DETECTED. 

IMPLICIT  REAL*8(A-H,L-Z,6) 

PRINT  1 , IFIN 

FORMATC/1X, 'FIN  CODE  IFIN=',I6/) 

INDUCE  ERROR  IN  ORDER  TO  GENERATE  TRACING  OF  CALLING  SUBROUTINES. 
X=-l . DO 
Y=X+DSQRTCX) 

IFCIFIN . LE. 0)  GO  TO  100 
STOP 
RETURN 


JuSrOuTInE 


MARCH 

SUBROUTINE  NUMBER  5 

IMPLICIT  REAL#3( A-H, L~Z»  $ ) 

COMMON  /VECS/XFC 101 ) , RMFC 101 ) , RPFC101 ) , MFC  101 ) ,MUF(101 ) , 

1  TETAF(101),BF(101), 

2  XN( 101 ),RMN(101),RPN(101),MN(101), MUN( 101), 

3  TETANdOl).BN(lOl),  XTEMP(  101 ) 

COMMON/THICKY/XTH( 1 002), TH( 1002) 

COMMON  /GAMA/G , G1 , G2 , G3 , G4 , G5 , G6 , G7 , G8 , 69 , G1 0 , G1 1 , G12 , G1 3 , G1 4 , G1 5 
G16,G17,G18,G19,G20 


JET0453 

JET0454 

JET0455 

JET0456 

JET0457 

JET0458 

JET0459 

JET0460 

JET0461 

JET0462 

JET0463 

JET0464 

JET0465 


COMMON  /PAR/PAI , PAI2 , DEG, XC , YC, MEXIT , MFIN, YMAX, DYO , DY, DYNEXT, 


STAB, DELTA, PSI1, PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO , 
TETSYM, TETL IM, DDY , DYMAX 
COMMON  /STAG/RHOO,NO,PO,TO, A0,MD0T1 
COMMON  /IPAR/ JMAX,  ICFO  ,  ITERO  ,  KF,  KN ,  IM,  IP,  J, 

KF2. 1  DEL , JDEL , JYXI , JXI , ILEAD , ILEADF, KCLEAD 
COMMON  /ROW/YF , YN, DXF , DXN 

COMMON  /CHARAC/XCHARF( 92 ) , YCHARFC  92 ) , XCHARNC  92) , YCHARNC  92) , 
RMCARF( 92 ) , RPCARFC92) , RMCARN(92) , RPCARN(92) , 
TCHARF(92),TCHARN(92), MUCARF ( 92 ) , MUCARN( 92) , 
CSIGNNC  92 ) , CSIGNF( 92 ) , MCHARNC  92) , MCHARFC  92) , 
MCHARI ( 92 ) 


C 

C 


COMMON  / ICHARA/ KCHARP , KCHARM , KCHARO 


ADVANCE  FLOW  FIELD  FROM  YF  TO  YN 
IM=KF 
IP  =  KF 
YN=YF+DY 
KN=KFO 

SEMI-INVERSE  INTEGRATION  FOR  FAN  POINTS. 

CALL  SEMINV 

NEW  GRID  POINTS  (JUST  INVERSE  MARCHING). 

CALL  GRIDN 

LOAD  FLOW  VARIABLES  FROM  SEMI-INVERSE  INTEGRATION  INTO  VECTORS 
CALL  LOADC 

CHARACTERISTIC  SCHEME  INTEGRATION  FOR  INNER  POINTS  (INVERSE  MARCH). 
CALL  INVMAR 
RETURN 
END 
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SUBROUTINE  INVMAR 
SUBROUTINE  NUMBER  6 

IMPLICIT  REAL*8CA-H, L-Z, $) 

COMMON  /VECS/XFC 101 ) , RMF( 101 ) , RPFC 101 ) ,MFC 101 ) ,MUFC 101 ) , 

1  TETAFC101),BFC101), 

2  XNC 101 ) , RMNC 101 ) , RPNC101),MNC101),  MUN(101 ) , 

3  TETANC101),BNC101), XTEMP( 101 ) 

COMMON/THICKY/XTHC 1002), THC 1002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7 , G8,G9,G10,G11,G12,G13»G14,G15 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PA I , PAI2, DEG, XC, YC, MEXIT , MFIN, YMAX, DYO, DY, DYNEXT , 

1  STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO, 

2  TETSYM,TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1  KF2, I DEL , JDEL , JYXI , JXI , ILEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF,YN, DXF, DXN 

INTEGRATION  WITH  INVERSE  CHARACTERISTICS  FOR  NEW  P0INTCX4, Y4) . 

OLD  POINTS  ARE  CXI , Y1 ) , (X2, Y2) . 

XI  IS  OBTAINED  BY  INVERSE  C-  FROM  X4 
X2  IS  OBTAINED  BY  INVERSE  C+  FROM  X4 
NOTE  THAT  XI  MAY  BE  NEGATIVE  CE.  G.  WHEN  X4=0). 

KN1=ILEAD-1 

IFCKN1 .LE.O)  CALL  FINC601) 

DO  1000  1=1, KN1 
14  =  1 

X4*XN( I ) 

Y4=YN 

IF4=CIM+IP)/2 

CALL  INTERPC  0 , IF4 , KF, X4, XF, RM4, RMF, RP4,  RPF) 

CALL  RFUNC(RM4,RP4,M4,MU4, TETA4) 

M14=M4 

MU14=MU4 

TETA14=TETA4 

M24=M4 

MU24=MU4 

TETA24=TETA4 

Y1=YF 

Y2=YF 

Y14=(Y1+Y4)/2.D0 
Y24=(Y2+Y4)/2. DO 
Xl=l .DIO 
X2=l .DIO 
RM4=1.D10 
RP4=1 .DIO 
ITER=0 
GO  TO  2 

CORRECTOR 

ITER=ITER+1 

AVERAGED  PROPERTIES  ON  C-( 14 ) , C+C 24 )  CHARACTERISTICS. 
RM14=(RMl+RM4)/2 . DO 
RP14=(RP1+RP4)/2.D0 
RM24= ( RM2+RM4 )/2 . DO 
RP24=(RP2+RP4)/2.D0 

M14,MU14,TETA14,  M24 , MU24 , TETA24  AVERAGED  ON  C-,C+  CHARACTERISTICS. 
CALL  RFUNCC RM14, RP14,M14,MU14,TETA14) 

CALL  RFUNCC RM24 , RP24 , M24 , MU24 , TETA24 ) 

CONTINUE 
NEW  XI, X2 
X10=X1 
X20=X2 

X1=X4-DY/DTANCTETA14-MU14) 

X2  =  X4-DY/DTANC  TETA24+MU24 ) 

IFCX2.LT. 0. )  CALL  FINC670) 

D14=DSQRTCCX1-X4)**2+DY**2) 

D24  =  DSQRTC  CX2-X4)**2+DY**2) 

INTERPOLATE  OLD  DISTRIBUTION  FOR  RM1,RP1,  RM2,RP2  AT  XI, X2. 

CALL  INTERPC0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 

CALL  INTERPC 0, IP , KF , X2 , XF , RM2 , RMF , RP2 , RPF) 
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C  NO  NEED  FOR  RE-AVERAGING  SINCE  IT  INTRODUCES  ONLY  HIGHER  ORDER 
C  CHANGES  INTO  THE  ITERATION  SCHEME. 

C  INTEGRATE  THE  CHARACTERISTIC  EQUATIONS  FOR  RM4,RP4  AT  X4,Y4. 

RM40SRM4 

RP40*RP4 

RM4=RMl+DELTA*DSIii(TETA14)*D14/<M14XY14) 

RP4=RP2+DELTAXDSIN(TETA24)XD24/(M24XY24) 

C  CONVERGENCE  TEST 

EPS=CDABSCX1-X10)+DABS<X2-X20))/DY+DABSCRM4-RM40)+DABS(RP4-RP40) 

IF(ITER.GT.ITERO)  GO  TO  10 

IFC  EPS .GT . EPSIL )  GO  TO  1 

RMN( I ) =  RM4 

RPNC I ) =  RP4 

CALL  RFUNCIRM4, RP4, MN< I ) , MUN< I ) , TETANC I ) ) 

1000  CONTINUE 
RETURN 

10  CONTINUE 

PRINT  11,I4,KN,IF4,IM,IP,KF,ITER,ITER0,EPS,EPSIL,X1,X2,X4,M14,M24 

11  FORMAT! IX, ' SUBR.  INVMAR .  14, KN, IF4, IM, IP, KF, ITER, ITER0= • ,815/ 

1  IX, *EPS,EPSIL,X1,X2,X4,M14,M24=',7D14.6/) 

CALL  FIN(611 ) 

RETURN 

- 5UTO3UTrNE~STMrNV - 

C  SUBROUTINE  NUMBER  7 

IMPLICIT  REALx8( A-H, L-Z, $ ) 

COMMON  /VECS/XF( 101 ) , RMF( 101),RPF(101 ) ,MF( 101 ) ,MUF( 101 ) , 

1  TETAF(101),BF(1Q1)> 

2  XN( 101 ) , RMN( 101 ) , RPN( 101 ) ,MN( 101) ,MUN( 101 ) , 

3  TETAN(101),BN(101), XTEMP( 101 ) 

COMMON/THICKY/XTH! 1002) ,TH(1002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2, DEG,XC,YC,MEX1T ,MFIN, YMAX, DYO, DY , DYNEXT , 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG, EPSIL, NUO, 

2  TETSYM.TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1  KF2, IDEL , JOEL , JYXI , JXI , ILEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF , YN , DXF , DXN 

COMMON  /CHARAC/XCHARF( 92 ) , YCHARF( 92) , XCHARNC 92) , YCHARNC  92) , 

1  RMCARF( 92) , RPCARFC  92 ) , RMCARNC  92) , RPCARNC  92) , 

2  TCHARF(92),TCHARN(92)>  MUCARFC  92 ) , MUCARN!  92) , 

3  CSIGNNC  92) , CSIGNF( 92) ,MCHARN(92) , MCHARFC92) , 

4  MCHARIC92) 

COMMON  /ICHARA/KCHARP , KCHARM, KCHARO 


JET0577 
JET0578 
JET0579 
JET0580 
JET 0581 
JET0582 
JET0583 
JET0584 
JET0585 
JET0586 
JET0587 
JET0588 
JET0589 
JET0590 
JET0591 
JET0592 
JET0593 
JET0594 
JET0595 
JET0596 
JET0597 
JET0598 
JET0599 
JetO^oo 

JET0601 
JET0602 
JET0603 
JET0604 
JET0605 
JET0606 
JET0607 
JET0608 
JET0609 
JET0610 
JET0611 
JET0612 
JET0613 
JET 06 14 
JET0615 
JET0616 
JET0617 
JET0618 
JET0619 
JET0620 
JET 06 21 
JET0622 
JET0623 
JET0624 
JET0625 
JET0626 
JET0627 
JET0628 
JET0629 
JET0630 
JET0631 
JET0632 
JET0633 
JET0634 
JET0635 
J  ET  06  36 
J  ET  06  37 
JET0638 
J  ET 06  39 
JET0640 
JET0641 
JET0642 
JET0643 
JET0644 
JET0645 
JET0646 


•-  A 
'A 


& 

& 


rv* 

& 


£ 

-V 


c 

c 

c 

c 


c 

c 

c 


COMPUTE  NEW  POINT  <X4,Y4),  BY  PASSING  A  C+  CHARACTERISTIC 
THROUGH  OLD  POINT  (X2.Y2).  BOTH  POINTS  ARE  ON  CHARACTERISTIC  LINE 
NUMBER  KC. 

IM=1 

DO  100  KC=1, KCHARO 
IF(CSIGNNCKC) . EQ . 0 . )  GO  TO  100 

PREDICTOR 

Y1=YF 

Y2=YF 

Y4=YN 

Y14=(Y1+Y4)/2.D0 

Y24=(Y2+Y4)/2.D0 

X2=XCHARF(KC) 

RM2=RMCARF(KC) 

RP2=RPCARF( KC ) 

M2=MCHARF(KC) 

MU2=MUCARF(KC) 

TETA2=TCHARF(KC) 

M14=M2 

MU14=MU2 

TETA14=TETA2 


FILE.  JETPR 
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TETA24STETA2 
X4*l . DIO 
Xl=l . DIO 
RM4S1 . DIO 
RP4=1 .DIO 
ITER=0 
GO  TO  2 
C 

C  CORRECTOR 
C 

1  ITER=ITER+1 

C  AVERAGED  PROPERTIES  ON  C-(14),C+<24)  CHARACTERISTICS. 
RM14=(RMl+RM4)/2. DO 
RP14=(RP1+RP4)/2.D0 
RM24=(RM2+RM4)/2.DO 
RP24=(RP2+RP4)/2.D0 

C  M14,MU14,TETA14,  M24,MU24, TETA24  AVERAGED  ON  C-,C+  CHARACTERISTICS. 
CALL  RFUNCCRM14, RP14,M14,MU14, TETA14) 

CALL  RFUNCCRM24,RP24,M24,MU24,TETA24) 

2  CONTINUE 
C  NEW  X4,X1 

X40=X4 

XIO=Xl 

X4-X2+DY/DTAN(TETA24+CSIGNF(KC)*MU24) 

X1=X4-DY/DTAN(TETA14-CSIGNF(KC)XMU14) 

D14=DSQRT((X1-X4)XX2+DYX*2) 

D24=DSQRT((X2-X4)XX2+DYXX2) 

C  INTERPOLATE  OLD  DISTRIBUTION  FOR  RM1.RP1,  AT  XI. 

CALL  INTERP(0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 

IF(J.GT.l)  GO  TO  22 
IFtCSIGNF(KC) . LT . 0 . )  GO  TO  22 
RP1=RP2 
22  CONTINUE 

C  NO  NEED  FOR  RE-AVERAGING  SINCE  IT  INTRODUCES  ONLY  HIGHER  ORDER 
C  CHANGES  INTO  THE  ITERATION  SCHEME. 

C  INTEGRATE  THE  CHARACTERISTIC  EQUATIONS  FOR  RM4.RP4  AT  X4,Y4. 
RM40=RM4 
RP40=RP4 

IFCCSIGNF(KC) .LT.O. )  GO  TO  21 

RM4=RM1+DELTA*DSIN(TETA14)*D14/CM14*Y14> 

RP4=RP2+DELTA*DSINCTETA24)*D24/(M24KY24) 

GO  TO  20 
21  CONTINUE 

RM4=RM2+DELTA*DSIN(TETA24)*D24/(M24XY24) 

RP4=RP1+DELTAXDSIN(TETA14)XD14/(M14XY14) 

20  CONTINUE 
C  CONVERGENCE  TEST 

EPS=( DABS <X4-X40)+DABS( XI -X10) )/DY+ DABS (RM4-RM40)+DABS(RP4-RP40) 
IF(ITER.GT.ITERO)  GO  TO  10 
IF( EPS .GT . EPSIL )  GO  TO  1 
CSIGNN(KC)=CSIGNF(KC) 

IFCX4.GT.0.J  GO  TO  30 
RMSAVE=RM4 
RM4=RP4+TETSYM 
RP4=RM4-TETSYM 
CSIGNN( KC) =-l . DO 
30  CONTINUE 

RMCARN(KC)=RM4 
RPCARNC  KC ) =RP4 

CALL  RFUNC(RM4,RP4,M4,MU4,TETA4) 

TCHARN( KC) =TETA4 
XCHARNC  KC) =DABS( X4 ) 

YCHARN( KC ) =Y4 
MUCARN( KC) =MU4 
MCHARN( KC) =M4 
100  CONTINUE 
RETURN 

10  CONTINUE 

PRINT  1 1, KC, KCHARO, IM, KF, ITER, ITERO, EPS, EPSIL, XI, X2,X4,M14,M24 

11  FORMAT ( IX, ' SUBR .  SEMINV.  KC, KCHARO , IM, KF, ITER, ITERO= ', 615/ 

1  IX, 'EPS, EPSIL, XI, X2,X4,M14, M24= ' , 7  D14 . 6/ ) 

CALL  FIN(711 ) 
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[NE  RFUNCCRM,RP,M,MU,TETA) 

SUBROUTINE  NUMBER  8 

IMPLICIT  REALXSCA-H,L-Z,8) 

COMMON  /VECS/XFC101),RMFC101),RPFC101),MFC101),MUFC101), 

1  TETAFC101),BFC101), 

2  XNC 101 ) , RMN( 101 ) # RPN( 101), MN (101 ) , MUNC 1 01 ) , 

3  TETANC101),BNC101),XTEMPC10J) 

COMMON/THICKY/XTHC 1002) ,THC1002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2, DEG, XC, YC, MEXIT, MFIN, YMAX, DYO, DY, DYNEXT , 

1  STAB, DELTA, PSI 1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL ,  NUO , 

2  TETSYM,TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO ,NO,PO,TO,AO, MDOT1 
COMMON  /IPAR/JMAX, KFO,  ITERO,KF,KN,  IM,  IP,  J, 

1  KF2, IDEL , JOEL , JYXI, JXI , ILEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF, DXN 
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10 


COMPUTE  M, MU, TETA  AT  A  POINT,  AS  FUNCTION  OF  RIEMANN  INVAR.  RM,RP. 
TETA*CRM-RP)/2. DO+TETLIM 
NU  =CRM+RP)/2 . DO 

NU=NU0-(G5*ARCTAN(G5XQ)-ARCTAN(Q)),  WHERE  Q=CMXX2-1 )**C-l/2) 

FIND  Q(NU) ,  AND  HENCE  M(NU),  THROUGH  NEWTON  RAPHSON  ITERATIONS. 
Q=-CNU-NU0VCG4-1 .DO) 

IFCQ.LE.O.)  CALL  FIN(801) 

ITER=0 

ITER=ITER+1 

QF=Q 

DNUDT=-CG4-1.DO)/C(1.DO+G4XQKX2)X(1 .D0+QXX2)) 
DNU=NU-(NU0-(G5XDATAN(G5XQ)-DATAN(Q) ) ) 

Q=Q+DNU/DNUDT 
IFCQ.LE.O.)  CALL  FINC811) 

EPS=DABSCQ-QF)/Q 
IF(  ITER  .GT .  HERO )  GO  TO  10 
IF( EPS .GT . EPSILXl . D-3 )  GO  TO  1 
M=DSQRT ( 1 . DO+1 . DO/QXX2) 

MU=DARSINC 1 . DO/M) 

RETURN 
CONTINUE 
CALL  FINC810  ) 

RETURN  llfTCR  p 


END 


Subroutine  INTERPc jneW, i , i<grid, x, xvec, rm, rmvec, rp, rpvec) 

SUBROUTINE  NUMBER  9 

IMPLICIT  REAL*8(A-H,L-Z,$) 

COMMON  /VECS/XF(101),RMF(101),RPF<101),MFC101),MUFC101), 

1  TETAF(101),BF(101), 

2  XN( 1 01 ) , RMN( 1 01 ),RPN(101),MN(101), MUNC 101), 

3  TETAN(101),BN(101), XTEMPC 101) 

COMMON/THICKY/XTHC1 002 ),THC 1002) 

COMMON  /GAMA/G, G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI, PAI2, DEG, XC.YC.MEXIT, MFIN, YMAX, DYO, DY, DYNEXT, 

STAB, DEL TA, PSI 1, PSIF, ZETA1, SIGMA, FRACG, EPSIL, NUO, 
TETSYM,TETLIM, DDY, DYMAX 
/STAG/RHOO, NO, PO, TO, AO, MDOT1 
/IPAR/JMAX,  KFO,  HERO ,  KF,  KN ,  IM,  IP,  J, 

L  KF2,IDEL, JDEL, JYXI, JXI, ILEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF, DXN 
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COMMON 

COMMON 


C 

C 

C 
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DIMENSION  XVECC 1 ) , RMVECC 1 ) , RPVECC 1 ) 


FIND  I  SUCH  THAT  XVECC I ). L E . X . AND . XVECC 1+1 ). GE . X 
FIND  RM, RP  BY  LINEAR  INTERPOLATION. 

NOTE  THAT  X  MAY  BE  NEGATIVE. 

IFCDABSCX) .LE. XVECC  KGRI D) )  GO  TO  901 
PRINT  900 , X, KGRI D, XVECC  KGRI D) 

F0RMATC/1X, D15.7,I10,4X,D15.7/) 

CALL  FINC  900 ) 

901  CONTINUE 

KG2=2XKGRID 
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FILE*  JETPR 

FORTRAN  A1 

I0*MIN0d,  KGRID-2) 

JET0793 

IC0UNT*0 

JET0794 

1 

I»IO 

JET0795 

SIGN1*1 . DO 

JET0796 

IF(I.GE.l)  GO  TO  10 

JET0797 

SIGN1*-1 . DO 

JET0798 

I*-I+2 

JET0799 

10 

CONTINUE 

JET0800 

IFd.GT.KGRID)  CALL  FIN(901) 

JET0801 

XXl*SIGNl*XVECd ) 

JET0802 

11*1 

JET 08 03 

IFCXX1 . LE . X)  GO  TO  11 

JET0804 

10*10-1 

JET 08 05 

IC0UNT*IC0UNT+1 

JET0806 

IF( ICOUNT . GT . KG2)  CALL  FIN(911) 

JET0807 

GO  TO  1 

JET0808 

11 

CONTINUE 

JET0809 

1*10+1 

JET0810 

SXGN2*1 . DO 

JET0811 

IFd.GE.l)  GO  TO  12 

JET0812 

SIGN2*-1 . DO 

JET0813 

I=-I+2 

JET0814 

12 

CONTINUE 

JET0815 

IFd.GT.KGRID)  CALL  FIN<912) 

JET0816 

XX2=SIGN2*XVECd) 

JET0817 

12*1 

JET0818 

IFCXX2.GE.X)  GO  TO  13 

JET0819 

10=10+1 

JET0820 

IC0UNT*IC0UNT+1 

JET0821 

IFUC0UNT.GT.KG2)  CALL  FINL913) 

JET0822 

GO  TO  1 

JET0823 

13 

CONTINUE 

JET0824 

F1=(XX2-X)/(XX2-XX1) 

JET0825 

j 

F2*l . D0-F1 

JET0826 

J 

IFCF1.LT.0.)  CALL  FINC991) 

JET0827 

IF( F2 . LT . 0 . )  CALL  FINC992) 

JET0828 

J 

RM1=RMF(I1) 

JET0829 

RP1=RPF( 11 ) 

J  ET 0830 

RM2=RMF( 12 ) 

JET0831 

RP2=RPF(I2) 

JET0832 

* 

IFCSIGN1 . LT . 0 . )  RM1=RPF(I1 l+TETSYM 

JET0833 

» 

IFLSIGNl . LT . 0 . )  RP1=RMF( 11 )-TETSYM 

JET0834 

■ 

IFCSIGN2 . LT . 0 . )  RM2=RPF( 12 l+TETSYM 

J  ET 0835 

, 

IFCSIGN2 . LT . 0 . )  RP2=RMFCI2)-TETSYM 

JET0836 

RM=F1*RM1+F2*RM2 

J  ET 0837 

RP=F1*RP1+F2*RP2 

JET0838 

RETURN 

END 

ItJTERK 
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bUBKUU  I  INt  INI  EKXl  JNEW,  1 1  ,  VAKU ,  VAK, KUK1 u, XU , XVfC) 

SUBROUTINE  NUMBER  10 

IMPLICIT  REAL*8(A-H,L-Z,$) 

COMMON  /VECS/XF( 101 ) , RMF( 101), RPFdOl ) , MFdOl )  ,MUFd01 ) , 

1  TETAFd01),BFd01), 

2  XNd01),RMNd01),RPNd01),MNd01),MUNd01), 

3  TETANd01),BNd01),XTEMPd01) 
COMMON/THICKY/XTHU002),THU002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,GS,G9,G10,G11,G12,G13,G14,G15 
G16,G17,G18,G19,G20 


COMMON  /PAR/PAI ,PAI2, DEG, XC, YC,MEXIT,MFIN, YMAX, DYO, DY, DYNEXT, 


1 

2 


C 

C 


STAB,  DELTA,  PSI1,  PSIF,  ZETA1,  SIGMA,  FRACG, EPSIL , NUO> 
TETSYM, TETL IM, DDY , DYMAX 
COMMON  /STAG/RHOO, NO, PO, TO, AO, MDOT1 
COMMON  /IPAR/ JMAX, KFO,ITERO,KF,KN,IM, IP, J, 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF,KCLEAD 

COMMON  /ROW/YF,YN, DXF, DXN 
DIMENSION  VAR(1),XVEC(11 

FIND  XO  AND  II  SUCH  THAT  XVECdl )<XO<XVEC( I 1+1 ) ,  AND  XO  CORRESPONDS 
TO  THE  LOCATION  AT  WHICH  VARO  IS  A  LINEAR  INTERPOLATION  OF  VAR(I). 
X0=1 .D23 
IFIRST=I1 

IF(Il.GT.O)  GO  TO  10 
I  FIRST  =  KGRID-IABS( 1 1 )  +  2 
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JETPR 
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10 

CONTINUE 

JET0865 

DO  1  II=IFIRST , KGRID 

JET0866 

I  =  II 

JET0867 

IF(Il.GT.O)  GO  TO  11 

JET0868 

I=KGRID-II+2 

JET0869 

11 

CONTINUE 

JET0870 

IF(I.LE.O)  CALL  FIN(lOOl) 

JET 08 71 

IFCI.GT. KGRID)  CALL  FINC1Q02) 

JET0872 

IF(I.EQ.l)  GO  TO  1 

JET0873 

IF(CVAR(I)-VARO)*(VARCI-l)-VARO).GT.O.)  GO  TO  1 

JET0874 

IFCVAR(I) .EQ.VAR(I-l))  GO  TO  1 

JET0875 

F1=(VAR( I )-VAR0)/( VARC I )-VAR( 1-1 ) ) 

JET0876 

F2=l . D0-F1 

JET0877 

IFCF1.LT.0.)  CALL  FIN(lOll) 

JET0878 

IFCF2.LT.0.)  CALL  FIN(1012) 

JET0879 

X0=F1XXVEC(I-1)+F2XXVEC(I) 

JET0880 

11=1-1 

JET0881 

GO  TO  2 

JET0882 

1 

CONTINUE 

JET0883 

2 

CONTINUE 

JET0884 

RETURN 

END 

JET0885 

JET0886 

3U7KPUTTNE  BKhAH 
SUBROUTINE  NUMBER  11 

IMPLICIT  REAL*8CA-H,L-Z,$) 

REAL  MB, MX, MY 

COMMON  /VECS/XF(101),RMF(101),RPF<101),MF<101),MUFC101), 

1  TETAF(101),BF(101), 

2  XN( 101 ) , RMNC101 ) , RPN( 101 ) , MN( 101 ) ,MUN( 101 ) , 

3  TET AN 1101), BN (101), XTEMPC 101) 

COMMON/THICKY/XTHU 002),  TH(  1002) 

COMMON  /GAMA/6, G1 , G2, G3 , G4 , G5 , G6 , G7 , G8 , G9 , G10, 611, G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PA I , PAI2, DEG, XC, YC, MEXIT ,MFIN, YMAX, DYO, DY, DYNEXT , 

1  STAB, DELTA, PSI1 , PSIF, 2ETA1 , SIGMA, FRACG, EPSIL ,  NUO , 

2  TETSYM, TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO ,NO,PO,TO,AO, MD0T1 
COMMON  /IPAR/JMAX, KFO ,ITERO,KF,KN,IM,IP,J, 

1  KF2, I DEL , JDEL , JYXI , JXI , ILEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF, DXN 


i 


i 


K 

I 


c 

c 


STORE  IN  BN( I ) . 


i 


COMPUTE  THE  BREAKDOWN  PARAMETER  AT  (1-1/2, K-l/2) . 
YB=0.5D0*(YF+YN) 

DYY = DY 
IM=2 

DO  1  1=2, KN 
X1=XN(I-1) 

X2=XN( I ) 

DXX=X2-X1 

IF(X2.GT.XF(KF) )  GO  TO  2 

CALL  INTERP(0,IM,KF,X1,XF,RM1,RMF,RP1,RPF) 

CALL  INTERPC  0,IM,KF,X2,XF, RM2 , RMF, RP2 , RPF ) 

CALL  RFUNC(RM1,RP1,M1,MU1,TETA1) 

CALL  RFUNCC RM2, RP2, M2, MU2 , TETA2 ) 
MX=0.5D0*C(MN(I)-MN(I-1))+(M2-M1) )/DXX 
MY=0 . 5D0XC (MN( I ) -M2 ) +( MN( I -1 )-Ml ) )/DYY 
MB=0.25D0*(MN(I-1)+MN(I)+M1+M2) 
TETAB=0.25D0*(TETAN(I-1)+TETAN(I)+TETA1+TETA2) 
G0REM=MBX*2*(1 . D0+G1*MB**2)**(G6-1 .DO) 

GRAD=MXXDCOS( TETAB)+MY*DSIN( TETAB ) 

B=G20*GQREMXGRAD 
GO  TO  3 
B=1 . D22 
BMC  I  )  =  B 
CONTINUE 
BN( 1 ) =BN(2 ) 

RETURN 
END 

— smwuuiiNt  UPAUY - 

SUBROUTINE  NUMBER  12 

IMPLICIT  REAL*8(A-H,L-Z,$) 

COMMON  /VECS/XFC 101),RMF(101),RPF(101),MF(101),MUF(101), 


TEITJ88  /" 
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JET0891 
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JET0918 
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JET0920 
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JET0922 
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JET0924 
JET0925 
JET0926 
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JET0928 
JET0929 
JET0930 
JET0931 
JET0932 
TE 1  0933 
JET0934 
JET0935 
JET0936 


FILE >  JETPR 
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C 

C 


1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3  TETAN(101),BN(101), XTEMP( 101) 

COMMON/THICKY/XTH( 1002) >TH(1002) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14, 
1  G16 , G17 , G18, G19, G20 

COMMON  /PAR/PAI , PAI2, DEG, XC, YC, MEXIT, MFIN, YMAX, DYO, DY, DYNEXT, 

1  STAB, DELTA, PSI1 , PSIF, ZETA1, SIGMA, FRACG, EPSIL , NUO, 

2  T  ETSYM , T  ET  L I M , DDY , DYMAX 
COMMON  /STAG/RHOO,NO,PO,TO,AO,MDOT1 
COMMON  /IPAR/JMAX, KFO, ITERO, KF, KN, 1M, IP , J, 

1  KF2,  IDEL ,  JDEL,  JYXI,  JXI,  HEAD,  ILEADF,  KCLEAD 

COMMON  /ROW/YF,YN, DXF, DXN 


1 

11 


COMPUTE  THE  MOLECULAR  THICKNESS  AT  END  POINTS  OF  EACH  ROW. 
IM=2 

XTH(J)*XF(KF) 

TH( J  )  =  0 . 

DTH0=N0*SIGMA*DY 
IFCJ.EQ.l)  GO  TO  11 
J1=J-1 

DO  1  JJ=1,J1 
XXI -XTH( J J ) 

CALL  INTERP( 0 , IM, KF, XXI , XF, RM1 , RMF, RP1, RPF) 

CALL  RFUNCCRM1 , RP1 ,M1 ,MU1 , TETA1 ) 

G0REM=1.D0+G1*M1**2 
DTH=DTH0/G0REM**G6 
THCJJ)=TH(JJ)+DTH 
CONTINUE 
CONTINUE 

RETURN  PLUMES 


END 


bUVRUU'f  l'NE""PrUMET> - 

SUBROUTINE  NUMBER  13 

IMPLICIT  REAL*8(A-H, L-Z, $) 

REALK4  XPL  YPL 

COMMON  /PLUME/XPL (1002,10) » YPL ( 1002) 

COMMON  /IPLUME/KPL, ITYPL (10) 

DIMENSION  VPLC92) 

COMMON  /VECS/XF< 101 ) , RMF( 101 ) , RPFC 101 ) , MFC101 ) , MUFC 101 ) , 

1  TETAF(101),BF(101), 

2  XNC101), RMN< 101),RPN(101),MN(101),MUN(101), 

3  TETAN(1Q1),BN(101), XTEMPC 101 ) 

COMMON/TH I CKY/ XT H< 1002),TH(1002) 

REALX4  YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON  /THICKX/YXIC20),XIC101,20) , XIPMC1QI , 20 ) , XIGRPC101 , 20 ) 

1  ,XIAPP(101,20),XIF(101,20) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2 , DEG, XC , YC, MEXIT , MFIN, YMAX, DYO , DY, DYNEXT, 

1  STAB, DELTA, PSI 1 , PSI F, ZETA1 , SIGMA, FRACG, EPSIL, NUO, 

2  TETSYM,TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO , NO , PO , TO , AO , MD0T1 
COMMON  /IPAR/JMAX, KFO , ITERO , KF, KN, IM, IP , J , 

1  KF2,IDEL, JDEL, JYXI, JXI, I  LEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF,YN, DXF, DXN 

COMMON  /CHARAC/XCHARFC  92 ) ,YCHARF(92), XCHARNC  92 ) , YCHARNC  92 ) , 

1  RMCARFI 92 ) , RPCARFI 92 ) , RMCARNI 92 ) , RPCARNC  92 ) , 

2  TCHARFC92), TCHARNC  92 ) , MUCARFC92 ) , MUCARNC  92 ) , 

3  CSIGNN(92),CSIGNF(92) ,MCHARN(92) ,MCHARF(92) , 

4  MCHARI (92) 

COMMON  /ICHARA/KCHARP , KCHARM, KCHARO 
COMPUTE  SPECIAL  POINTS  AT  Y=YN,  AND  STORE  THEM  AS 
(XPLCJ, IPL) >YPL(J)=YN) . 

J  IS  THE  MARCHING  INDEX  OF  YN. 

IPL  =  1 , 2 , . . . , KPL  IS  THE  "PLUME"  INDEX.  PRESENTLY  KPL.LE.5 
VPLCIPL)  IS  A  VALUE  DEFINING  THE  "PLUME"  CURVE. 

ITYPL ( IPL )  IS  THE  TYPE  OF  CURVE.  IT  DEFINES  CURVES  AS  FOLLOWS: 

DO  NOTHING 

REAL  PLUME.  IT  IS  THE  BREAKDOWN  SURFACE,  DEFINED 
BY  A  CONSTANT  VALUE  OF  THE  BREAKDOWN  PARAMETER  B. 
SET  VPL ( IPL ) =B . 


JET0937 
JET0938 
JET0939 
JET0940 
G15, JET0941 
JET0942 
JET0943 
JET0944 
JET0945 
JET0946 
JET0947 
JET0948 
JET0949 
JET0950 
JET0951 
JET0952 
JET0953 
JET0954 
JET0955 
JET0956 
JET0957 
JET0958 
JET0959 
JET0960 
JET0961 
JET0962 
JET0963 
JET0964 
JET0965 
JET0966 
JET0967 
JET0968 


- JET096T' 

JET0970 
JET0971 
JET0972 
JET 097 3 
JET0974 
JET0975 
JET0976 
JET0977 
JET0978 
JET0979 
JET0980 
JET0981 
JET0982 
JET0983 
G15, JET0984 
JET0985 
JET0986 
JET0987 
JET0988 
J  ET  0989 
JET0990 
JET0991 
JET0992 
JET0993 
JET0994 
JET0995 
JET  0996 
JET0997 
JET0998 
JET0999 
JET1000 
JET1001 
JETI002 
JET 100 3 
JETI004 
JET1005 
JET1006 
JET1007 
JET1008 


ITYPL ( IPL ) =0 
ITYPL ( IPL ) =1 
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Pi 


W 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ITYPL( IPL )*2 
ITYPL(IPL)S3 


ITYPL(IPL)a4 


ITYPL (IPL )a5 


CONSTANT  MACH-NUMBER  LINE.  VPL(IPL)aM. 

A  SINGLE  STREAMLINE.  VPL(IPL)  IS  SET  TO  THE  EXIT 
X-COORDINATE  OF  THAT  STREAMLINE. 

A  SINGLE  C+  CHARACTERISTIC  LINE  STARTING  AT  THE  CORNER 
VPL (IPL)  IS  SET  TO  THE  INDEX  KC  OF  THAT  CHARACTERISTIC 
LINE. 

A  CONSTANT  LATERAL  (X)  OPACITY  LINE.  VPL(IPL)  IS  SET 
TO  THE  VALUE  OF  THE  (CONSTANT)  OPACITY. 


2001 


DEFINE  ITYPL(IPL)  AND  VPL(IPL) 

KPL=10 

IF(KPL.GT.IO)  CALL  FIN<1301) 

DO  2000  IPLal,KPL 

GO  TO  (2001,2002, 2003,2004,2005, 2006, 2007,2008, 2009, 2010),  IPL 
ITYPL(IPL)a4 


2002 


2003 


2004 


2005 


2006 


2007 


2008 


2009 


2010 


2000 


C 

C 

C 

C 


200 


VPL(IPL)al 
GO  TO  2000 
ITYPL(IPL)a4 
VPLdPL  )=KCHARP 
GO  TO  2000 
ITYPL ( IPL ) =4 
VPL(IPL)=19 
GO  TO  2000 
ITYPL(IPL)=4 
VPL( IPL )=31 
GO  TO  2000 
ITYPL(IPL)=4 
VPL(IPL)a47 
GO  TO  2000 
ITYPL(IPL)=4 
VPL(IPL )=55 
GO  TO  2000 
ITYPL ( IPL ) =1 
VPL (IPL )=0 . 02D0 
GO  TO  2000 
ITYPLdPL  )  =  1 
VPL ( IPL ) =0 . 03D0 
GO  TO  2000 
ITYPL ( IPL ) =1 
VPL ( IPL ) =0 . 05D0 
GO  TO  2000 
ITYPL ( IPL ) =1 
VPL(IPL)=0.08D0 
GO  TO  2000 
CONTINUE 

COMPUTE  "PLUME"  POINTS  AT  Y=YN 
DO  1000  IPL  =1 , KPL 
ITYP=ITYPL ( IPL ) 

IF(ITYP.EQ.O)  GO  TO  1000 
GO  TO  (1,2, 3, 4, 5),  ITYP 
CONTINUE 

BREAKDOWN  SURFACE  PLUME. 

NOTE  THAT  DUE  TO  DIFFERENCE-CENTERING  OF 
Y-COORDINATE  IS  0.5*(YF+YN),  RATHER  THAN 
IN  THE  PLOTTING  CODE. 

BO=VPL ( IPL ) 

XTEMP( 1 )  =XN( 1 ) 

DO  11  1=2, KN 

XTEMP( I)=0.5DQ#(XN(I)+XN(I-1)) 

CONTINUE 
1  =  2 

CALL  INTERXd,  I ,  BO ,  BN,  KN,  XBO ,  XTEMP ) 

XPL ( J , IPL ) =XBO 
GO  TO  1001 
CONTINUE 

FIND  BY  INTERPOLATION  THE  X-COORDINATE  WHERE  M=MPL 
IF(J.GT.l)  GO  TO  200 
XPL (J,IPL)=XC 
GO  TO  1001 
CONTINUE 
MPL  =VPL ( IPL ) 

I  =  -KN 


GRADIENTS,  THE 
YN.  IT  CAN  BE 


ACCURATE 

ADJUSTED 


34 


JET1009 

JET1010 

JET1011 

JET1012 

JET1013 

JET1014 

JET1015 

JET1016 

JET1017 

JET1018 
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JET1021 

JET1022 

JET1023 

JET1024 

JET1025 

JET1026 

JET 1027 

JET1028 

JET1029 

JET1030 

JET1031 

JET1032 

JET1033 

JET1034 

JET1035 

JET1036 

JET1037 

JET1038 

JET1039 

JET1040 

JET1041 

JET1042 

JET1043 

JET1044 

JET1045 

JET1046 

JET1047 

JET1048 

JET1049 

JET 1050 
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JET1052 
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JET1061 

JET1062 
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JET1066 

JET1067 

JET1068 
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JET1078 

JET  1 079 
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CALL  INTERXC1,I,MPL,MN,KN,XM0,XN) 

XPLC J, IPL)*XMO 
60  TO  1001 

3  CONTINUE 

:  STREAMLINE  INTERPOLATION. 

IF(J.GT.l)  GO  TO  300 
XPLCJ,IPL)*VPLCIPL) 

GO  TO  1001 

300  CONTINUE 
XSF»XPLCJ-1,IPL) 

ISF=2 

I5N*2 

CALL  INTERPC  0, ISF, KF,XSF, XF, RMSF, RMF, RPSF, RPF) 
CALL  RFUNCCRMSF,RPSF,MSF,MUSF,TETASF) 
XSN=XSF+DY*DTANCPAI2-TETASF) 

ITER=1 

301  ITER*ITER+1 

CALL  INTERP( 1 , ISN, KN, XSN, XN, RMSN, RMN, RPSN, RPN) 
CALL  RFUNC( RMSN, RPSN, MSN, MUSN,TETASN) 

TETAAV=0 . 5D0*CTETASF+TETASN) 
XSN=XSF+DY*DTANCPAI2-TETAAV) 

IF(ITER. LT. ITERO+2)  GO  TO  301 

XPL(J,IPL)=XSN 

GO  TO  1001 

4  CONTINUE 

:  CHARACTERISTIC  LINE. 

KC*IDINT(VPL(IPL)+1 .D-5) 

IF(J.GT.l)  GO  TO  41 
XPLC J» IPL)=XCHARF(KC) 

GO  TO  1001 
41  CONTINUE 

XPLCJ,IPL)=XCHARN(KC) 

IFCCSIGNNCKC) . EQ . 0 . )  XPLC J, IPL)*1 . E33 
GO  TO  1001 

5  CONTINUE 

:  CONSTANT  LATERAL  CX)  OPACITY 
CALL  OPACX 
XIC=VPLCIPL) 

DO  51  11=2, KF 

I1=KF-II+1 

12=11+1 

XI1=XICI1,JXI) 

XI2=XICI2,JXI) 

IFC  CXIC-XI1 )*CXIC-XI2) . GT . 0 . )  GO  TO  51 
F2=CXI2-XIC)/CXI2-XI1) 

Fl=l . D0-F2 

IFCF1.LT.0.)  CALL  FINC1351) 

IFCF2.LT.0.)  CALL  FINC1352) 

XIFC=F2#XFC I1)+F1*XFC 12) 

GO  TO  52 

51  CONTINUE 
XIFC=1 . D30 

52  CONTINUE 
XPLCJ,IPL)=XIFC 
GO  TO  1001 

1001  CONTINUE 
IFCJ.GT.l)  GO  TO  1002 
YPL ( J ) =YC 

GO  TO  1000 

1002  CONTINUE 
YPL  C  J ) =YN 

1000  CONTINUE 
RETURN 
END 
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SUBROUTINE  NUMBER  14 

IMPLICIT  REAL*8C A-H, L-Z, ♦ ) 

REALX4  XPL  YPL 

COMMON  /PLUME/XPL C 1002, 10), YPL C 1002) 

COMMON  /IPLUME/KPL, ITYPLC10) 

COMMON  /VECS/XFC 101 ) , RMFC 101), RPFC 1 01 ) ,MFC 101 ) , MUFC 101 ) , 
1  TETAFC101),BF(101), 


JET1081 
JET1082 
JET1083 
JET1084 
JET1085 
JET1086 
JET1087 
JET1088 
JET1089 
JET1090 
JET1091 
JET1092 
JET1093 
JET1094 
JET1095 
JET1096 
JET1097 
JET1098 
JET1099 
JET1100 
JET1101 
JET1102 
JET1103 
JET1104 
JET1105 
JET1106 
JET1107 
JET 11 08 
JET1109 
JET1110 
JET1111 
JET1112 
JET1113 
JET1114 
JET1115 
JET1116 
JET1117 
JET1118 
JET1119 
JET1120 
JET1121 
JET 11 22 
JET1123 
JET1124 
JET1125 
JET1126 
JET1127 
JET1128 
JET1129 
JET1130 
JET1131 
JET1132 
JET1133 
JET1134 
JET1135 
JET1136 
JET1137 
JET1138 
JET1139 
JET1140 
JET1141 
JET1142 
JET1143 
JET1144 


J u I ilij 

JET1146 
JET 11 47 
JET1148 
JET1149 
JET1150 
JET1151 
JET1152 
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FORTRAN  A1 


JETPR 


2  XNd01),RMN(101),RPN(101),MNd01),MUN(101),  JET1153 

3  TETANdOl )  ?  BN(  101 ) »XTEMP(101 )  JET1154 
COMMON/THICKY/XTH(1002),TH(1002)  JET1155 
COMMON  /GAMA/G, 01, G2,G3,G4,G5, G6,G7,G8, G9,G10, Gil, 012, G13, G14,G15, JET1156 


1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT, 

1  STAB,  DELTA,  PSI1 , PSIF, ZETA1 ,  SIGMA,  FRACG, EPSIL,NU0, 

2  TETSYM,TETLIM,DDY,DYMAX 
COMMON  /STAG/RHOO , NO , PO , TO , AO , MDOT1 
COMMON  /IPAR/JMAX, KFO, ITERO, KF, KN, IM, IP, J, 

1  KF2,IDEL,  JOEL,  JYXI,JXI,  HEAD,  ILEADF, KCLEAD 

COMMON  /ROW/YF,YN, DXF, DXN 

COMMON  /CHARAC/XCHARFC 92 ) , YCHARF( 92) , XCHARNC 92) , YCHARNt 92 ) , 

1  RMCARFC  92) , RPCARF( 92) , RMCARN( 92) , RPCARN<  92) , 

2  TCHARFC  92) ,TCHARN( 92) , MUCARF( 92) , MUCARNC92) , 

3  CSIGNNC92) ,CSIGNF(92) ,MCHARN(92),MCHARF(92) , 

4  MCHARK92) 

COMMON  /ICHARA/KCHARP, KCHARM, KCHARO 
DIVIDE  LINE  Y=YN  INTO  KN-1  INTERVALS. 

THE  X-GRID  IS  NON-UNIFORMLY  DEFINED  AS  FOLLOWS! 

(1)  (XCHARN( I ) , YCHARN( I ) ) ,  (XCHARFd ) , YCHARFC I ) ) ,  1=1,2, .. . , KCHARP, 
DENOTE  NEW  AND  OLD  (FORMER)  CHARACTERISTIC  <C+)  POINTS.  LET  1=1 
AND  I=KCHARP  CORRESPOND  TO  THE  LEADING  AND  BOUNDARY 
CHARACTERISTICS  CC+). 

THE  GRID  CONSISTS  OF  TWO  SEGMENTS.  THE  SO-CALLED  FLAT  SEGMENT 
IS  BETWEEN  X=0  AND  X=XLEAD=XCHARN(KCLEAD) .  THE  SECOND  IS  THE 
FAN  SEGMENT.  IT  IS  FROM  XLEAD  TO  XBOUND=XCHARN(KCHARP) . 


(2) 
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(3)  THE  FAN  SEGMENT  IS  INITIALLY  DIVIDED  INTO  FRACGX(KFO-l)  INTERVAL SJET1 180 


(4) 


DEFINED  BY  THE  FAMILY  OF  C+  CHARACTERISTIC  LINES  MCHARI(l)  TO 
MCHARI (KCHARP). 

THE  FLAT  SEGMENT  IS  DIVIDED  INTO  (1-FRACG)*(KF0-1)  EQUAL 
INTERVALS,  AS  LONG  AS  THEY  ARE  NOT  SMALLER  THAN  THE  AVERAGE 


WHEN  THEY  ARE,  THEIR  NUMBER  IS  REDUCED,  BUT  NOT 


(5) 


1 

11 


FAN  INTERVAL 
BELOW  THREE. 

KCLEAD  IS  INITIALLY  1.  IT  IS 
IS  AT  LEAST  TWICE  THE  AVERAGE 
ILEADF=ILEAD 
KCLEAD=0 

DO  1  KC=1, KCHARP 
IF(CSIGNN(KC) . LT . 0 . )  GO  TO  1 
KCLEAD=KC 

KFAN=KCHARP-KCLEAD 
XL EAD=XCHARN( KCLEAD) 

XBOUND=XCHARN( KCHARP) 
DX1=(XB0UND-XLEAD)/DFL0AT(KFAN) 

I F(XL EAD/DX1 . GT . 2 . DO )  GO  TO  11 

CONTINUE 

CONTINUE 

IF( KCLEAD . EQ .  0)  CALL  FIN(1401) 

IF(KCLEAD.EQ.  KCHARP)  CALL  FINU402) 
ILEAD=IDINT(XLEAD/DXl)+2 
IFdLEAD+KFAN.GT.KFO)  ILEAD  =  KFO-KFAN 
ILEAD1 =ILEAD-1 
KN=ILEAD+KFAN 

IF(KN.GT.KFO)  CALL  FIN(1411) 
DX=XLEAD/DFLOAT ( ILEAD1 ) 

XN( 1 )  =  0  . 

DO  2  1=1 , ILEAD1 

XN( I ) =XN(1 )+DX#DFLOAT( 1-1 ) 

CONTINUE 

DO  3  I=ILEAD, KN 

XN( I )=XCHARN(KCLEAD+I-ILEAD) 

CONTINUE 
RETURN 
END 


UPDATED  SO  THAT  THE 
FAN  INTERVAL. 


FLAT  SEGMENT 


ysr£P 


Subroutine  ystep 

C  SUBROUTINE  NUMBER  15 

IMPLICIT  REAL#8(A-H, L- 
REALX4  XPL , YPL 

COMMON  /PLUME/XPL( 1002, 10), YPL (1002) 

COMMON  /IPLUME/KPL , ITYPL( 10) 

COMMON  /VECS/XF( 101 ) , RMF( 101 ) ,RPF(101),MF(101),MUF(101), 
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"JFITZTr 
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FORTRAN  A1 


C 

C 

C 

C 


1  TETAF(101),BF(101),  JET1225 

2  XN(101),RMNC101),RPN(101),MN(101),MUNC101),  JET1226 

3  TETAN(101),BNfl01), XTEMPC 101 )  JET1227 

COMMON/THICKY/XTHC 1002) *  TH( 1002)  JET1228 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, JET1229 

1  G16,G17,G18,G19,G20  JET1230 

COMMON  /PAR/PAI,PAI2, DEG, XC,YC,MEXIT,MFIN,YMAX,DYO,DY, DYNEXT,  JET1231 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA,  FRACG, EPSIL,NUO,  JET1232 

2  TETSYM,TETLIM,DDY, DYMAX  JET1233 

COMMON  /STAG/RHOO, NO , PO, TO, A0,MD0T1  JET1234 

COMMON  /IPAR/JMAX, KFO, ITERO,KF,KN,IM,IP,J,  JET1235 

1  KF2, 1  DEL ,  JOEL ,  JYX1 ,  JXI ,  HEAD,  ILEADF,KCLEAD  JET1236 

COMMON  /ROW/YF,YN, DXF, DXN  JET1237 

COMMON  /CHARAC/XCHARFC92) ,YCHARF(92) , XCHARNC92) ,YCHARN(92) ,  JET1238 

1  RMCARFC92) , RPCARFC  92) , RMCARN( 92) , RPCARN(92) ,  JET1239 

2  TCHARFC  92) , TCHARNC  92) ,MUCARF( 92) , MUCARNC  92) ,  JET1240 

3  CSIGNNC  92) , CSIGNFC92) , MCHARN( 92) , MCHARFC92) ,  JET1241 

4  MCHARIC92)  JET1242 

COMMON  /ICHARA/KCHARP , KCHARM, KCHARO  JET1243 

COMPUTE  NEXT  Y-STEP.  JET1244 

DYNEXT  IS  DEFINED  AS  THE  MINIMAL  "TRIANGULATION"  Y-STEP  DYT ,  0BTAINEDJET1245 
BY  FORWARD  INTERSECTION  OF  C-,C+  CHARACTERISTICS  FROM  ADJACENT  GRID  JET1246 


POINTS  XI, X2. 

DYMIN=1.D40 
DO  1  1=3, KF 
X1=XFCI-1) 

X2=XF( I ) 

DX=X2-X1 

TP1=DTAN(TETAF( 1-1 )-MUF( 1-1 ) ) 

TP2=DTAN(TETAF( I )+MUF( I ) ) 

F1=-TP2/(TP1-TP2) 

IFCFl.LE.O.)  PRINT  555, I , XI , X2,DX,TP1 ,TP2, FI 
555  F0RMATC/1X,  '  I ,  XI ,  X2,  DX,  TP1 ,  TP2,  Fl  =  M5, 6D14 . 6/) 
IFCFl.LT.O.)  CALL  FINC1501) 

DYT=F1*DX*TP1 

IF( DYT . LE. 0 . )  CALL  FINC1502) 

DYMIN=DMIN1 C  DYMIN, STABXDYT ) 

1  CONTINUE 

DYNEXT=DYMIN 

renetdurn  mve 


SUBROUTINE  MOVE - 

SUBROUTINE  NUMBER  16 

IMPLICIT  REAL*8(A-H,L-Z,$) 

REALX4  XPL , YPL 

COMMON  /PLUME/XPLC 1002, 10), YPL (1002) 

COMMON  /IPLUME/KPL , ITYPl ( 10 ) 

COMMON  /VECS/XFC 101),RMF(101),RPF(101),MF(101),MUF(101), 
TETAFClOl)rBF(lOl), 

XN( 1 0 1 ) , RMNC 101),RPN(101),MN(101),MUN(101), 
TETAN(101),BN(101), XTEMP( 101 ) 

COMMON/THICKY/XTHC 1002), TH(1 002) 

COMMON  /GAMA/G, G1 ,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14, 
G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2, DEG , XC , YC, MEXIT , MFIN, YMAX, DYO, DY, DYNEXT, 
STAB, DELTA, PSI 1 , PSI F, ZETA1 , SIGMA, FRACG, EPSIL,NUO, 
TETSYM, TETLIM, DDY, DYMAX 
COMMON  /STAG/RHOO, NO, PO, TO, AO, MD0T1 
COMMON  /IPAR/JMAX, KFO , ITERO , KF , KN , IM, IP, J, 

1  KF2,IDEL,JDEL,JYXI,JXI,ILEAD,ILEADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF, DXN 

COMMON  /CHARAC/XCHARFC  92 ) , YCHARFC  92 ) , XCHARNC 92) , YCHARN( 92) , 

1  RMCARFC92) , RPCARFC 92) ,RMCARN( 92), RPCARNC 92), 

2  TCHARFC92) ,TCHARN( 92), MUCARFC 92), MUCARNC 92), 

3  CSIGNN(92),CSIGNF(92), MCHARNC  92) ,MCHARF(92), 

4  MCHARI ( 92 ) 

COMMON  /ICHARA/KCHARP , KCHARM, KCHARO 
STORE  NEW  LINE  (N)  IN  OLD  LINE  (F). 

KF  =  KN 

KF2=2XKF 

YF=YN 

DO  1  1=1, KN 
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XF<I)*XN(I) 

RMF(I)SRMN(I) 

RPFCI)=RPN(I) 

MF(I)*MN(I) 

MUF(I)SMUN(I) 

TETAFCI)=TETAN(I) 

BF(I)=BN(I) 

CONTINUE 
DO  2  KC=1 , KCHARO 
IFCCSIGNNCKC) . EQ . 0 . )  GO  TO  2 
XCHARFC KC) =XCHARNC  KC) 
YCHARFCKC)=YCHARNCKC) 

RMCARF( KC) =RMCARN( KC) 
RPCARFCKC)=RPCARNCKC) 
TCHARFCKC) =TCHARNCKC) 
MUCARF(KC)=MUCARNCKC) 
MCHARFCKC)=MCHARNCKC) 

CSIGNFC  KC)=CSIGNNCKC) 

CONTINUE 

RETURN 

END  _ 
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SUBROUTINE  NUMBER  17 

IMPLICIT  REAL*8(A-H,L-Z,$) 

REAL*4  XPL,YPL 

COMMON  /PLUME/XPL<1002,10),YPL<1002) 

COMMON  /IPLUME/KPL, ITYPLCIO) 

COMMON  /VECS/XF( 101 ) , RMF< 101 ) * RPFC 1Q1),MF( 101 ) ,MUFC101) , 

1  TETAFC101),BFC101), 

2  XNC 101) , RMN( 101 ) , RPN( 101 ) , MNC 101 ) , MUNC 101 ) , 

3  TETANa01),BNC101),XTEMP<101) 

COMMON/THICKY/XTH( 1002), TH( 1002) 

REAL*4  YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON  /THICKX/YXI(20),XIC101,20),XIPM(101,20),XIGRP(101,20) 

1  , XI APPC1 01, 20 ),XIF(101,20) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2, DEG, XC, YC, MEXIT ,MFIN, YMAX, DYO , DY, DYNEXT , 

1  STAB, DELTA, PSI1 , PSIF, ZETA1, SIGMA, FRACG, EPSIL,NU0, 

2  TETSYM, TETLIM, DDY , DYMAX 
COMMON  /STAG/RHOO , NO , PO , TO , AO , MDOT1 

COMMON  /CHARAC/ XCHARFC 92 ) , YCHARFC 92 ) , XCHARNC 92 ) , YCHARNC 92 ) , 

1  RMCARFC  92) , RPCARFC  92 ) , RMCARNC  92 ) , RPCARNC  92 ) , 

2  TCHARFC  92) , TCHARNC  92 ) , MUCARFC  92 ) , MUCARNC  92 ) , 

3  CSIGNN(92),CSIGNF(92) ,MCHARN(92) ,MCHARF(92) , 

4  MCHARIC92) 

COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP, J, 

1  KF2, 1 DEL , JDEL , JYXI , JXI , I  LEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF, DXN 
COMPUTE  X-OPACITY. 

BEGIN  FROM  LIMITING  CHARACTERISTIC  OF  AN  ASSUMED  P.M.  FAN. 

XIO  --  THE  THICKNESS  BETWEEN  THE  LIMITING  CHARACTERISTIC  AND  THE 
BOUNDARY  CHARACTERISTIC  OF  THE  NUMERICAL  COMPUTATION. 

DO  12  1=1, KFO 
XIFC I , JXI )=XF( I ) 

XI(I,JXI)=0. 

XIPMC I , JXI ) =0 . 

XIGRPC I , JXI )=0 . 

XIAPPC I , JXI ) =0 . 

I  CONTINUE 

IFCJ.EQ.l)  GO  TO  1000 
PSILIM=TETLIM 

XLIM=XC+( YF-YC )/DTAN(PSILIM) 

XBOUND=XFCKF) 

KPM=10 

DX=(XLIM-XBOUND)/DFLOATCKPM) 

SUM=0 . 

DO  1  1=1, KPM 

Xl=XBOUND+DFLOATCI-l)*DX 

X2=X1+DX 

PS1=PAI2-DATAN(  CX1-XO/C  YF-YC) ) 

PS2=PAI2-DATANC  CX2-XO/CYF-YC) ) 
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FORTRAN  A1 


Q1=(PS1-PSILIM)/G5 
Q2=(PS2-PSILIM)/G5 
IF(I.EQ.KPM)  Q2*l . D-IO 
IF(Q2.LT.O. )  CALL  FIN(1701) 

F1=G11X( DSIN(Q1 ) )XX(2 . DO/(G-I .  DO ) ) 

F2=G11X(DSIN(Q2) )xx<2 . DO/(G-l . DO) ) 

SUM=SUM+DXX( Fl+F2)/2 . DO 
CONTINUE 

XIO*SUMX(NOXSIGMA) 

RE-EVALUATE  XIO  FOR  A  RING-JET. 

IF(DELTA.EQ.O.)  GO  TO  14 
M=MFIN 

CALL  MFUNC(M,F,ETA,TETA) 

PSI*TETA+DARSIN(1 . DO/M) 

GOREM=l . DO+G1XMXX2 
GOR=MXX2-1.DO 
CALL  HINTERCM, HM) 

DELTOB»0.5DOXDSQRT(GOR)X(1.DO/(MEXITXETA)+DSIN(TETA)/M)/DSIN(PSI) 
1  +G15XHM/2.D0 

EVER-SIGMAXNOXYC/(MXDSIN(TETA)XDSIN(PSI)XGOREMXXG6) 

GGG=2. DO-DEL TO BX(G+1 . DO )/2 . DO 
IFCDABS(GGG) .GT.l .D-IO)  GO  TO  15 
PRINT  16,  DELTOB,G,GGG 

>  F0RMAT(/1X, 'FROM  OPACX.  GGG  NEARLY  ZERO.  EXPRESSION  FOR  XIO  IS', 

1  IX, 'SINGULAR.  DELTOB, G, GGGS ',3012.4/) 

CALL  FIN( 1715) 

5  CONTINUE 

EVER3 EVER/ GGG 

XI 0 =EVERX ( ( YF/YC ) XXGGG-1 . DO )/ ( YF/ YC ) 

»  CONTINUE 

XICKF, JXI)=XIO 

XIPM(KF, JXI )=XIO 

XIGRP(KF, JXI )=XIO 

KF1=KF-1 

DO  2  11*1, KF1 

I=KF-II+1 

X1=XF(1) 

X2=XF(I-1) 

DX=X1-X2 

F1=1.D0/(1.D0+G1XMF<I  )XX2)XXG6 

F2=l . D0/( 1 . D0+G1XMF( 1-1 )XX2)XXG6 

DTNUM=(N0XSIGMA)XDXX(F1+F2)/2.D0 

XI(I-1,JXI)=XI(I , JXI )+DTNUM 

XIPM(I-1,JXI)=1.D24 

XIGRP(I-1, JXI)=1 .D24 

PS1=PAI2-DATAN((X1-XC)/(YF-YC)) 

PS2=PAI2-DATAN((X2-XC)/(YF-YC)) 

IFCPS2.GT.PSI1)  GO  TO  2 

Ql=(PSl-PSILIM)/05 

Q2=(PS2-PSILIM)/G5 

IF(Q1 .  LT.  0  . )  CALL  FINU711) 

F1=G11X(DSIN(01) )xx(2.D0/(G-l.D0)) 
F2=G11X(DSIN(Q2))XX(2.D0/(G-1.D0)) 

DTPM=(NOXSIGMA)XDXX(F1+F2)/2.DO 
XIPMCI-1,  JXI)=XIPM(I,  JXD  +  DTPM 
DISTl=DSQRT((Xl-XC)xx2+(YF-YC)xx2) 
DIST2=DSQRT((X2-XC)x*2+(YF-YC)xx2) 

XC1 =KCLEAD+I -I L  EAD 
KC2=KC1-1 

IF(KC2 . LT . KCLEAD)  GO  TO  21 
Ml =MCHARI( KC1 ) 

M2=MCHARI(KC2) 

CALL  MATCHd  ,  Ml  ,MG1 ,  MOBI 1  ,MABI1 ) 

CALL  MATCH( 1-1 , M2, MG2,M0BI2, MABI2 ) 

Fl=l .00/(1 .D0+G1XMG1XX2)X*G6 
F2=l . D0/(1 . D0+G1XMG2XX2)*XG6 
DTGRP=(N0XSIGMA)XDXX(F1+F2)/2.D0 
XIGRPCI-1,JXI)=XIGRP(I,JXI)+DTGRP 
L  CONTINUE 
CONTINUE 

APPROXIMATE  THICKNESS  XIAPPC I , JXI ) .  BASED  ON  CLOSED-FORM  INTEGRATION 
DO  3  1=1, KF 
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XIAPPCI, JXI)=1.D24 
KC=KCLEAD+C I-ILEAD) 

IFCDELTA . EQ . 0 . )  GO  TO  3 
IFCKC.LT. KCLEAD)  GO  TO  3 
IFCXFCI) .LT.XCHARFC1))  GO  TO  3 
M*MCHARICKC) 

CALL  MFUNCCM,F,ETA,TETA) 

PSI=TETA+DARSINCl.DO/M) 

GOREM=1.DO+G1XMXX2 
GOR=Mxx2-l.DO 
CALL  HINTERCM,HM) 

DELT0B=0.5D0KDSQRT(GOR)X(l.D0/<MEXIT*ETA)+DSINCTETA)/M)>'DSIN(PSI) 
1  +G15XHM/2.D0 

EVER=SIGMAXNOXYC/CMXDSINCTETA)XDSINCPSI)XGOREMXXG6) 

GGG=2 . DO-DELTOBXCG+1 . D0)/2 . DO 
IFC DABSCGGG) . GT . 1 . D-l 0 )  GO  TO  25 
PRINT  26,  I , KC, M, DELTOB, G, GGG 

FORMATC/1X, 'FROM  OPACX.  GGG  NEARLY  ZERO.  EXPRESSION  FOR  XIO  IS', 

1  IX, ‘SINGULAR.  I, KC,M=‘, 15,012.4/ 

2  IX, ‘DELTOB, G,GGG=',3D12. 4/) 

CALL  FINC 1725) 

CONTINUE 

EVER=EVER/GGG 

XIAPPCI ,  JXI  )  =  EVERXC  CYF/YOXXGGG-1 .  DO  )/CYF/YC) 

3  CONTINUE 
1000  CONTINUE 
RETURN 
END 
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- SUBROUTINE  LOATH? 

SUBROUTINE  NUMBER  18 

IMPLICIT  REALX8CA-H, L~Z, ♦) 

REALX4  XPL , YPL 

/PLUME/XPLC 1002, 10 ),YPLC 1002) 

/IPLUME/KPL , ITYPL  CIO) 

/VECS/XFC 101 ) , RMFC 1 01 ) ,RPFC101)»MFC101),MUFC1Q1)> 

L  TETAFCIOD.BFCIOI), 

l  XNC 1 01 ) , RMNC 101),RPNC101),MNC101),MUNC101), 

5  TETANC101),BNC101), XTEMPC 101 ) 

COMMON/THICKY/XTHC 1002 ) , THC 1002 ) 

REALX4  YXI , XI , XI PM, XI GRP ,XIAPP,XIF 

COMMON  /THICKX/YXI C20),XIC101,20),XIPMC101,20), XIGRPC 101,20) 
l  ,XIAPPC101,20),XIF(101,20) 

COMMON  /GAMA/G,G1,G2,G3.G4,G5.G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
l  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI,PAI2, DEG , XC, YC , MEXIT, MFIN , YMAX, DYO , DY, DYNEXT , 

STAB, DELTA, PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL , NUO, 
TETSYM, TETLIM, DDY, DYMAX 
/STAG/RHOO , NO, PO, TO, AO, MD0T1 
/IPAR/JMAX, KFO, ITERO,KF,KN,IM,IP,J> 

KF2, 1  DEL , JDEL , JYXI , JXI , IL  EAD, I LEADF, KCLEAD 
/ROW/YF,YN, DXF, DXN 

/CHARAC/XCHARFC  92) , YCHARFC  92 ) , XCHARNC  92) , YCHARNC92) , 
RMCARFC92) , RPCARFC  92) , RMCARNC  92 ) , RPCARNC  92) , 
TCHARFC92 ) ,TCHARNC92), MUCARFC  92) , MUCARNC  92) , 

CSIGNNC  92 ) , CSIGNFC  92 ) , MCHARNC  92) , MCHARFC92) , 

MCHARIC  92) 

/ICHARA/KCHARP, KCHARM, KCHARO 

VARIABLES  OF  GRID  POINTS  IN  THE  FAN  SEGMENT  FROM  THE 
SEMI-INVERSE  INTEGRATION  CIN  SUBR.  SEMINV).  NOTE  THAT  GRID  POINTS 
XNC I )  WERE  ALREADY  DETERMINED  IN  SUBR.  GRIDN. 

DO  1  I=ILEAD, KN 

KC=KCLEAD+I-IIEAD 

IFCKC . GT . KCHARP )  CALL  FINC1801) 

RMNC I )=RMCARNCKC) 

RPNC I )= RPCARNC  KC) 

MNCI) =MCHARNCKC) 

MUNC I ) =MUCARNC  KC) 

TETANC I ) =TCHARNC  KC) 

CONTINUE 

gg;uRN  hiuFunc 
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SUBROUTINE  NUMBER  19 

IMPLICIT  REAL*8£A-H,L-Z,$) 

REAL**  XPL,YPL 

COMMON  /PLUME/XPL(10Q2,10),YPL(1002) 

COMMON  /IPLUME/KPL , ITYPL ( 10) 

COMMON  /VECS/XF(101),RMF(101),RPF(101),MF(101),MUF(101), 

1  TETAF(101),BF(101), 

2  XN( 101 ) , RMN( 1 01 ) , RPN( 1 01 ) ,MN( 101) ,MUN( 101 ) , 

3  TETAN( 101 ) ,  BN( 101 ) , XTEMP( 101 ) 

COMMON/THICKY/XTH( 1002),  TH( 1002) 

REALX4  YXI,XI,XIPM,XIGRP,XIAPP,XIF 

COMMON  /THICKX/YXI (20),XIC101,20),XIPM(101,2Q), XIGRPC 101 > 20 ) 

1  , XIAPP( 101 , 20),XIF(101,20) 

COMMON  /GAMA/G, G1 , G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2, DEG, XC, YC, MEXIT, MFIN, YMAX, DYO, DY, DYNEXT, 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,  EPSIl,NU0, 

2  TETSYM,TETLIM, DDY , DYMAX 
COMMON  /STAG/RHOO ,N0,P0,T0,A0, MDOT1 
COMMON  /I PAR/ UMAX, KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2, 1 DEL , JDEL , JYXI , JXI , I  LEAD, ILEADF, KCLEAD 

COMMON  /ROW/YF, YN, DXF, DXN 

COMMON  /CHARAC/XCHARF(92) , YCHARFC 92) , XCHARNC92) , YCHARN(92) , 

1  RMCARFC92) , RPCARFC92) , RMCARN(92) , RPCARNC92) , 

2  TCHARF(92) , TCHARN( 92) , MUCARFC92) , MUCARN( 92) , 

3  CSIGNN( 92) , CSIGNFC  92) ,MCHARN(92) ,MCHARFC92) , 

4  MCHARI (92) 

COMMON  /ICHARA/KCHARP, KCHARM, KCHARO 
COMPUTE  NU  AS  FUNCTION  OF  MACH  NUMBER  M.  NOTE  THAT  THE  P.M. 
DEFINITION  OF  NU  HAS  BEEN  MODIFIED  BY  ADDING  A  CONSTANT.  THE  USUAL 
CHOICE  OF  THE  CONSTANT  IS  SUCH  THAT  NU=0  FOR  INFINITE  M. 

Q  =  1 .  D0/DSQRT(M**2-1 .  DO ) 

NUFUNC=NU0-(G5*DATAN(G5*Q)-DATAN(Q)) 

?nETdURN  MS£T 
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SUBROUTINE  HMSeT 
SUBROUTINE  NUMBER  20 

IMPLICIT  REAL*8(A-H,L-Z,$) 

REALX8  KAPAOB 

COMMON  /GAMA/G, G1.G2,G3,G4,G5,G6,G7,G8,G9,G10, Gil, G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PAI2, DEG, XC, YC, MEXIT , MFIN, YMAX, DYO , DY , DYNEXT , 

1  STAB, DELTA, PSI 1 , PSI F , ZETA1 , SIGMA, FRACG, EPSIL,NU0, 

2  TETSYM, TETLIM, DDY, DYMAX 
COMMON  /IPAR/JMAX,KFO,ITERO,KF,KN,IM,IP,J, 

1  KF2.IDEL, JDEL, JYXI, JXI, I  LEAD, ILEADF, KCLEAD 

COMMON  /GRP/DMINV,MHINV(101 ) ,HMV(101 ) 

COMMON  /IGRP/KHM 

ROUTINE  FOR  THE  C+  DERIVATIVE  DUE  TO  RING  SYMMETRY  (GRP). 

KHM=51 

IF(KHM.GT.lOl)  CALL  FIN(2001) 

MINVO=l .DO/MEXIT 
DMINV=MINVO/DFLOAT (KHM-1 ) 

M=MEXIT 
SUM=0. 

KHM1 =KHM-1 
DO  1  1=1, KHM1 
MF=M 

MHINV(I)=MINVO-DFLOAT(I-1)*DMINV 
M= 1 . DO/MHINV ( I ) 


DM=M-MF 
Ml =M-DM 
M2=M-DM/2.D0 
M3=M 

CALL  MFUNC(M1,F1,ETA1,TETA1) 

CALL  MFUNCCM2, F2,ETA2,TETA2) 

CALL  MFUNCCM3, F3,ETA3,TETA3) 

SUM=SUM+DMX(Fl+4.D0xF2+F3)/6 

ETA=ETA3 

TETA=TETA3 

PSI=TETA+DARSIN(1 .DO/M) 
NORM=((3.DO-G)/4 ,D0)X(MX*2-1 


DO 


DO )**0 . 75D0/ 


41 
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FORTRAN  A1 


PA 


ty, 


K 


11 


12 
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1  CDSIN(PSI)X(l.DO+GlXMXX2)XXG14) 

HM=SUMXNORM 

HMV(I)*HM 

GOREM-1.DO+G1KMXX2 

GOR=MX*2-1.DO 

DELTOB=0.5DOXDSQRT(GOR)X(1.DO/(MEXITXETA)+DSIN(TETA)/M)/DS1N(PSI) 
1  +<(G+1.D0)/(2.D0X(3.D0-G)))XHM 

EPSIOB=DELTOB/DSQRT<GOR)-DSIN(TETA)/(MXDSIN(PSI)) 

KAPAOBsl . DO 

IF(DABS(PAI2-TETA).GT.l.D-6) 

1KAPA0B=DTAN(TETA)XEPSI0B 

IAMDOB=EPSIOB-DELTOBXGOREM/(GORXOSQRT(GOR>) 

PRINT  11,I,M,HM,TETAXDEG,PSIXDEG 

FORMATC/1X,  •  I,M,HM,TETA,PSI=',I5,5D12.4> 

PRINT  12, DELTOB, EPSIOBxDEG, KAPAOB»DEG, LAMDOB*DEG 
FORMAT (  IX, * DELTOB, EPSIOB, KAPAOB, LAMD0B= • , SX, 5D12 .4) 

CONTINUE 
MHINV(KHM)=0 . 

HMV(KHM) S1 . DO 

BSURM  MFUNC 


- SUBROUTINE  W-UNC(M7F-,  ETA, TE1A) - 

SUBROUTINE  NUMBER  21 

IMPLICIT  REALX8<A-H,t-Z,$) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16 ,  G17 ,  G18 ,  G19,  G20 

COMMON  /PAR/PAI,PAI2, DEG, XC,YC, MEXIT, MFIN,YMAX,DYO,DY,DYNEXT, 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG,EPSIL,NUO, 

2  NUPT1 , TETLIM 

COMMON  /IPAR/JMAX, KFO, ITERO, KF,KN, IM, IP, J, 

1  KF2, I DEL , JDEL , JYXI , JXI , ILEAD, ILEADF, KCLEAD 

COMMON  /GRP/DMI NV , MHINV ( 1 0 1 ) , HMV < 1 0 1 ) 


V 


L- 


vr 

& 


11 

1 


NU=NUFUNCCM) 

TETA=NUFUNC(MEXIT)+PAI2-NU 
GOREM=l .D0+G1XMXX2 
G0R=MXX2-1 .DO 

F=(MXX21X(G0REMXXG13)XDSIN(TETA)/G0RXX1.25DQ 
GOREMl=l . D0+G1*MEXITX*2 
GOR1=MEXIT**2-1.DO 

ETA=( CGOREM/GOREM1 )XXG14)X( (GOR1/GOR)XXO . 25D0) 

RETURN 

END 

SUBROUTINE  HINTER(M,H) 

SUBROUTINE  NUMBER  22 

IMPLICIT  REAL*8(A-H, L-Z, $) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,GI7,G18,G19,G20 

COMMON  /PAR/PAI , PAI2 , DEG, XC, YC, MEXIT , MFIN, YMAX, DYO , DY, DYNEXT, 

1  STAB,  DELTA,  PSI1 , PSIF, ZETA1 , SIGMA, FRACG, EPSIL, NUO, 

2  TETSYM, TETLIM, DDY,DYMAX 
COMMON  /IPAR/ JMAX, KFO, ITERO, KF,KN,IM, IP, J, 

1  KF2.IDEL, JDEL .JYXI, JXI, ILEAD, ILEADF, KCLEAD 

COMMON  /GRP/DMINV, MHINVL 1 01 ) , HMV( 1 01 ) 

COMMON  /IGRP/KHM 
COMPUTE  H(M)  BY  INTERPOLATION 
MINV=1 .DO/M 

I=KHM-IDINT(MINV/DMINV-1 .D-9)-l 
IF( I . GE . 1 . AND . I . LT . KHM)  GO  TO  1 
PRINT  11, I, KHM. M, MEXIT 

FORMATC/ IX , • I , KHM, M, MEXIT  =,,2I5,2D14.6/) 

CALL  FINC2201) 

CONTINUE 

F1=(MINV-MHINV(I+1))/DMINV 
F2=l . D0-F1 

IF( FI . IT . -1 . D-9 )  CALL  FIN(2210) 

IFCF2.LT.-1 . D-9)  CALL  FIN(2211) 

H=F1#HMV( I )+F2#HMV( 1+1 ) 

RETURN 
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FILE*  JETPR  FORTRAN  A1 


IMPLICIT  REALX8CA-H,  L-Z,#) 

COMMON  /VECS/XFC 1 01 ) , RMF( 101), RPF( 1 01 ) , MF( 101 ) , MUF( 101), 

1  TETAF(101),BF(101), 

2  XN(101),RMN(101),RPN(101),MN(101),MUN(101), 

3  TETANC 101 ) *  BNC 101 ) , XTEMP( 101 ) 

COMMON  /ROW/YF, YN, DXF, DXN 

COMMON  /GAMA/G,G1 ,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/PAI , PA 12, DEG, XC, YC,MEXIT ,MFIN, YMAX, DYO, DY, DYNEXT , 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG, EPSIL,NUO, 

2  T  ET  SYM , T  ET  L I M , D DY , DYMAX 
COMMON  /IPAR/JMAX, KFO, ITERO, KF, KN, IM, IP, J , 

1  KF2, IDEL, JDEL, JYXI , JXI , ILEAD, ILEADF, KCLEAD 

COMMON  /GRP/ DMI NV , MH I NV ( 1 0 1 ) , HMV ( 1 0 1 ) 

COMMON  /IGRP/KHM 

C  COMPUTE  H(M)  AND  THE  ALFA-DERIVATIVES 
M=MOB 

CALL  MFUNC(M, F, ETA,TETA) 

PSI=TETA+DARSINC1 .DO/M) 

CALL  HINTERCM, HM) 

GOREM=1.DO+G1*M**2 

GOR=M*»2-1.DO 

DELTOB=0.5DOKDSQRT(GOR)*(1.DO/(MEXITXETA)+DSINCTETA)/M)/DSIN(PSI) 

1  +G15XHM/2.D0 

FOB=(G7*GOREM)*XG2/M 
FAB=F0BX(YF/YC)XXDELT0B 
CALL  AREAFC FAB, MAB) 

C  COMPUTE  MABI  FROM  THE  INVERSE  PROBLEM  SOLUTION 
COTAV=(XF(I)-XC)/(YF-YC) 

PSIO=PAI2-DATANCCOTAV) 

EVY=YFXDLOG(YF/YC)/(YF-YC)-1.DO 

PSIN=PSIO 

DO  1  ITER=1,50 

PSI=PSIN 

M=DSQRT( 1 . D0+G4/DTAN( (PSI-TETLIM)/G5)XX2) 

M=DMAX1 CM, MEXIT) 

CALL  HINTERCM, HM) 

CALL  MFUNCCM, F, ETA,TETA) 

GOREM=1.DO+G1XMXX2 

GOR=MXX2-1.DO 

DELT0B=0.5D0XDSQRT(G0R)X(1.D0/CMEXITXETA)+DSIN(TETA)/M)/DSINCPSI) 

1  +G15XHM/2.D0 

EPSIOB=DELTOB/DSQRT(GOR)-DSINCTETA)/(MXDSINCPSI)) 

LAMDOB=EPSIOB-DELTOBXGOREM/(GORXDSQRTCGOR)) 

COTN=COTAV+LAMDOBXEVY/DSIN(PSI)XX2 

PSIN=PAI2-DATANCC0TN) 

DPSI=PSIN-PSI 

IFCDABSCDPSD.LT. l.D-9)  GO  TO  11 

I  CONTINUE 

PRINT  12,I,ITER,PSI,PSIN,DP5I,M,XFCI),YF,XC,YC 
12  F0RMATC/1X, • I , ITER, PSI , PSIN, DPSI ,M, XFC I ) , YF, XC, YC= '// 

1  IX, 2I4,8D11 . 3/) 

CALL  FINC2301) 

II  CONTINUE 

C  USING  MOBI=M  AS  COMPUTED  FROM  THE  INVERSE  PROBLEM,  FIND  MABI. 

MOBI=M 

M=M0BI 

CALL  MFUNCCM, F, ETA, TETA) 

PSI=TETA+DARSINC1 .DO/M) 

CALL  HINTERCM, HM) 

G0REM=1 .D0+G1XMXX2 
G0R=MX*2-1 . DO 

DELT0B=0 .5DOXDSQRTCGOR)XC1.DO/CMEXITXETA)+DSINCTETA)/M)/DSINCPSI) 

1  +G15XHM/2.D0 

FOB=CG7XGOREM)XXG2/M 
FAB=FOBXCYF/YC)XXDELTOB 
CALL  AREAFC FAB , MABI ) 

^URN  A9EAF 


SUBROUTINE  AREAFC F, M ) 
SUBROUTINE  NUMBER  24 

IMPLICIT  REALX8CA-H,L-Z,$) 
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COMMON  /GAMA/G, G1,G2,G3,G4,G5,G6,G7, G8, G9,G10,G11, 612,013,014, 615, JET1729 

1  G16,G17,G18.-G19,G20  JET1730 

COMMON  /PAR/PAI,PAI2,DEG,XC,YC,MEXIT,MFIN,YMAX,DY0,DY,DYNEXT,  JET1731 

1  STAB, DELTA, PSI1,PSIF,ZETA1, SIGMA, FRACG, EPSIL , NUO,  JET1732 

2  TETSYM, TETLIM, DDY, DYMAX  JET1733 

COMMON  /IPAR/JMAX, KFO , ITERO, KF, KN, IM, IP, J,  JET1734 

1  KF2, 1 DEL, JOEL , JYXI , JXI, I  LEAD, ILEADF, KCLEAD  JET1735 

COMMON  /GRP/DMINV,MHINV( 101 ) , HMVC 101 )  JET1736 

COMMON  /IGRP/KHM  JET1737 

C  COMPUTE  MACH  NUMBER  M  FROM  AREA  RATIO  FUNCTION  F  JET1738 

C  F=((2/CG+1))*<1+<G-1)*M*X2))X*<(G+1)/(2X(G-1)))/M  JET1739 

C  INITIAL  GUESS  IS  MIN  JET1740 

El=(FXMEXIT)x*(l.D0/G2)/G7  JET1741 

E2=( El-1 . D0)/G1  JET1742 

E3=DMAXl<E2,MEXITxx2)  JET1743 

MIN=DSQRT(E3)  JET1744 

EMN=MIN  JET1745 

DO  1  1=1,100  JET1746 

EMO=EMN  JET1747 

G0REM=1.D0+G1XEM0X»2  JET1748 

GOR=EMOX*2-1.DO  JET1749 

F0=(G7XG0REM)XXG2/EM0  JET1750 

DF=FO-F  JET1751 

PRINT  123,1, EMO, EMN, FO, F, DF, GOR, GOREM  JET1752 

FORMATdX,  ' I , EMO,  EMN,  FO,  F, DF, GOR, GOREM=  *,I5,7D12.4)  JET1753 

DFDM=FOXGOR/( EMOXGOREM)  JET1754 

DMN=DF/DFDM  JET1755 

EMN=EMO-DMN  JET1756 

EPSEM=DABS(DMN/EMN)  JET1757 

IF( EPSEM . LT . 1 . D-10)  GO  TO  11  JET1758 

I  CONTINUE  JET1759 

CALL  FINC2401)  JET1760 

II  CONTINUE  JET1761 

M=EMN  JET1762 

RETURN  JET1763 

END  JET1764 
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